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PREFACE 


The  objective  of  this  research  project  is  to  understand  phenomena  of  seismic  wave  generation 
and  propagation  that  affect  regional  and  teleseismic  seismograms  used  for  nuclear  test  mon¬ 
itoring.  This  annual  technical  report  includes  two  studies.  The  first  is  a  theoretical  study 
of  the  seismic  radiation  from  explosions  detonated  in  non-spherical  cavities,  i.e.,  cylindrical 
tunnels  of  finite  length  embedded  in  a  homogeneous,  isotropic,  elastic  medium.  A  frequency 
domain  boundary  element /discrete  wavenumber  algorithm  is  applied  to  compute  seismic 
wavefields  from  nuclear  and  non-nuclear  explosions  located  at  various  positions  along  a  tun¬ 
nel  axis.  The  computed  radiation  patterns  display  strong  shear  wave  generation  and  show 
significant  differences  between  the  two  types  of  explosions.  In  the  second  study,  we  experi¬ 
mentally  and  numerically  investigate  the  scattering  of  an  acoustic  P  wave  incident  on  a  highly 
irregular,  random  acoustic-elastic  interface  to  determine  whether  enhanced  backscattering 
occurs.  The  experiments  involve  ultrasonic  waves  reflected  from  a  glass  surfaced  etched  to 
produce  a  highly  irregular  3-D  surface.  We  find  that  2-D  numerical  simulations  predict  the 
3-D  experimental  results  well  at  small  incident  angles.  Both  numerical  and  experimental 
results  strongly  support  the  presence  of  enhanced  backscattering.  The  report  on  this  latter 
study  is  a  preprint  of  a  paper  submitted  to  the  Journal  of  the  Acoustical  Society  of  America. 
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WAVEFIELDS  FROM  AN  EXPLOSION  ALONG  THE  AXIS 


OF  A  FINITE  TUNNEL 


Summary 

This  study  examines  the  source  radiation  from  explosions  detonated  in  a  cylindrical  cavity 
of  finite  length  (tunnel)  embedded  in  an  infinite  homogeneous  elastic  isotropic  background 
medium.  The  radiated  wavefields  are  modeled  using  a  frequency  domain  boundary  ele¬ 
ment/discrete  wavenumber  algorithm.  The  algorithm  employs  the  indirect  formulation  of 
boundary  integral  equations  for  fluid  and  elastic  media.  For  axially  symmetric  problems,  the 
explosion  source  is  modeled  as  fictitious  surface  sources  distributed  on  the  cavity  boundary. 
Upon  discretizing  the  cavity  boundary  into  elements  of  uniform  source  distribution,  and  then 
imposing  boundary  conditions  at  each  element,  a  system  of  equations  is  obtained  with  the 
fictitious  source  distribution  on  each  element  as  the  unknowns.  Elements  of  the  resulting 
coefficient  matrix  are  integrals  of  displacement  and  stress  Green’s  functions  over  boundary 
elements.  Once  the  equivalent  boundary  sources  are  determined,  wavefields  inside  and  out¬ 
side  the  cavity  are  easily  calculated.  We  apply  the  algorithm  to  study  two  specific  cases  of 
explosion  sources:  nuclear  and  non-nuclear.  A  nuclear  explosion  is  specified  by  assigning  a 
very  high  compressional  velocity  of  10  km/sec  inside  the  cavity,  representing  loading  of  the 
cavity  wall  by  a  shock  wave.  An  ordinary  pressure  wave  traveling  in  air  (330  m/sec)  is  used 
for  a  non-nuclear  explosion.  The  results  show  different  source  radiation  patterns  between 
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the  two  types  of  explosions,  especially  when  the  explosion  is  located  off-center  in  which  case 
the  non-nuclear  explosion  radiation  displays  strong  directivity  effects.  In  contrast,  radiation 
from  a  nuclear  explosion  does  not  depend  on  the  source  position  inside  the  tunnel.  Both 
types  of  explosions  radiate  significant  shear  wave  energy  outside  the  cavity. 

Introduction 

Seismic  radiation  from  a  nuclear  explosion  is  generally  more  complicated  than  that  of  a 
simple  isotropic  point  source.  This  has  been  shown  by  various  observations  (e.g.,  Day  et 
al.,  1983;  Wallace  et  al.,  1983,  1985;  Priestley  et  al.,  1990).  The  observations  from  nuclear 
explosions  may  be  explained  by  one  or  more  of  the  following  physical  mechanisms:  (1)  the 
nuclear  explosion  itself;  (2)  tectonic  release;  (3)  spall;  (4)  anisotropic  and  heterogeneous 
media  near  the  source;  and  (5)  asymmetry  of  the  source.  So  far  no  single  one  of  these 
mechanisms  has  proven  sufficient  to  explain  all  of  the  observed  data  (Masse,  1981;  Gupta 
and  Blandford,  1983;  Johnson,  1988). 

The  purpose  of  this  study  is  to  investigate  the  seismic  radiation  from  an  explosive  source 
placed  in  a  non-spherical  cavity.  The  effects  of  spherical  cavities  on  seismic  explosions  was 
studied  in  the  early  sixties  (e.g.,  Latter  et  al.,  1961).  Glenn  et  al.  (1985)  and  Rial  and  Moran 
(1986)  studied  the  coupling  mechanism  and  the  radiation  pattern  from  an  explosion  in  an 
ellipsoidal  cavity  in  an  unbounded  medium.  Helmberger  et  al.  (1991)  showed  that  correct 
modeling  of  near-field  seismograms  observed  after  nuclear  explosions  requires  distinct  source 
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characteristics  as  well  as  the  local  <  rustal  structure.  They  modeled  the  source  characteristics 
by  using  asymmetric  sources  inside  ellipsoidal  cavities.  Lately,  Zhao  and  Harkrider  (1992) 
showed  in  their  theoretical  investigation  that  asymmetry  of  the  source  region  has  a  more 
significant  effect  on  shorter  period  radiation  than  longer  period  radiation.  In  this  study,  we 
have  investigated  the  generation  of  compressional  and  shear  waves  from  explosions  placed 
along  the  axis  of  cylindrical  tunnels  of  various  diameters. 

Modeling  far-field  radiation  from  sources  placed  in  a  cavity  poses  difficulties  for  the  finite 
difference  method  (FDM)  and  the  finite  element  method  (FEM),  because  of  large  scale 
differences  between  the  cavity  surface  and  formation  extent.  The  ability  to  calculate  the 
far-field  wavefields  is  severely  restricted  by  the  computer  memory  space  needed  by  these 
methods.  The  accuracy  of  these  methods  is  also  hampered  by  grid  dispersion  and  inaccurate 
handling  of  the  fluid/formation  interface.  For  a  finite  cavity,  the  commonly  used  discrete 
wavenumber  method  based  on  the  vertical  wavenumber  representation  is  no  longer  applicable. 
To  overcome  these  difficulties,  an  indirect  boundary  element  method  (BEM)  combined  with 
a  discrete  wavenumber  method  based  on  the  Sommerfeld  integral  representation  is  used. 

The  BEM  was  first  established  through  the  direct  formulation  by  Jawson  (1963)  and 
Symm  (1963)  for  potential  theory,  Rizzo  (1967)  for  elastostatics,  and  Cruse  (1968)  for  elas- 
todynamics.  Applications  of  the  indirect  BEM  in  elastodynamics  only  appeared  recently  for 
seismic  wave  scattering  by  surface  topographies  (e.g.,  Wong,  1982;  Kawase,  1988;  Sanchez- 
Sesma  and  Campillo,  1991).  Along  a  similar  line,  work  based  on  boundary  integral  equation 
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and  discrete  vertical  wavenumber  formulations  (Campillo  and  Bouchon,  1985;  Bouchon  et 
al .,  1989)  have  emerged.  Bouchon  and  Schmitt  (1989)  applied  this  method  in  studying  full 
waveform  acoustic  logging  in  an  infinite  irregular  borehole.  However,  the  formulation  of 
this  method  was  limited  to  an  infinite  borehole  in  an  infinite  homogeneous  medium.  More 
recently,  Bouchon  (1993)  and  Dong  (1993)  developed  a  BEM  algorithm  to  model  downhole 
seismic  sources  in  layered  isotropic  and  transversely  isotropic  media.  Dong  (1993)  also  used 
BEM  to  model  sources  in  semi-infinite  boreholes. 

In  this  study,  we  extend  the  BEM  method  for  sources  in  semi-infinite  boreholes  to  model 
wave  radiation  from  an  explosion  source  in  a  finite  cavity  embedded  in  a  homogeneous 
medium. 

Method 

Indirect  Integral  Formulation 

If  a  volume  point  source  is  placed  inside  a  cavity,  the  total  displacement  potential  in  the 
cavity  fluid  is  the  sum  of  a  direct  potential  due  to  the  source  and  a  reflected  potential  due  to 
the  boundary.  In  the  case  of  steady  state  radiation  (or  in  the  frequency  domain)  the  reflected 
field  can  be  expressed  as  an  integral  of  a  fictitious  source  distribution  over  the  cavity  surface, 
with  a  Green’s  function  being  the  integrand.  Therefore,  the  displacement  potential  in  the 
fluid  is 

6 i(x)  =  6i  +  f  dS"tfi(x,x')0(x')  for  x  G  Vb  +  B  ,  (1) 

J  B 
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where  the  volume  of  cavity  fluid  and  the  cavity  surface  are  denoted  by  Vf,  and  B,  respectively. 
Subscript  1  stands  for  the  fluid  region  and  <j>i  is  the  incident  potential.  The  fictitious  source 
distribution  over  the  cavity  surface  is  denoted  by  t/>(x').  The  integral  kernel,  gx,  is  the  scalar 
Green’s  function  in  an  infinite  homogeneous  medium  which  is  of  the  well-known  form 

eifc/|x-x'| 

S,(x'X')  =  4?i^'  <2> 

where  kj  =  u/c  is  the  wavenumber  for  the  fluid. 

The  above  representation  is  not  only  physically  intuitive  but  also  mathematically  rig¬ 
orous.  In  fact,  this  representation  can  be  obtained  from  the  mathematical  formulation  of 
Huygen’s  principle  with  assistance  of  the  uniqueness  and  equivalence  principles.  That  is,  the 
influence  of  the  elastic  medium  on  the  wavefield  inside  the  cavity  is  equivalent  to  impress¬ 
ing  a  sheet  of  fictitious  sources  on  the  boundary  between  the  fluid  and  the  elastic  medium. 
Fictitious  source  distribution  density  is  the  unknown  function  to  be  determined. 

For  the  source-free  elastic  medium  outside  the  cavity,  the  displacement  field  can  be  ex¬ 
pressed  as 

U  2(x)=  /<fS'G(x,x')-¥(x')  for  xeK  +  fl,  (3) 

J  B 

where  'I'(x')  is  a  vector  fictitious  source  distribution  on  the  boundary.  G(x,  x')  is  the  dyadic 
Green’s  function  for  displacement  and  has  the  following  form 

G(x,  x')  =  [k2plg0(x,  x')  +  VV  [g0(x,  x')  -  ga(x,  x')]}  •  (4) 

Here,  ga  and  g0  are  scalar  Green's  functions  of  the  same  form  as  in  (2)  except  that  kj  is 
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changed  to  ka  for  the  dilatational  wave  and  to  kp  for  the  shear  wave.  Equation  (3)  says  that 
the  displacement  in  the  region  outside  the  cavity,  (K)>  results  from  vector  fictitious  sources 
distributed  along  the  boundary.  For  axially  symmetrical  problems  in  cylindrical  coordinates, 
the  vector  fictitious  source  can  be  decomposed  into  source  distributions  in  radial  and  vertical 
directions.  Therefore,  we  have  two  unknown  distribution  functions  to  determine  in  order 
to  calculate  the  displacement  field  in  the  elastic  medium.  Uniqueness  of  the  calculated 
displacement  is  guaranteed  by  the  equivalence  principle  and  more  fundamentally  by  the 
uniqueness  theorem. 

Implementation  of  the  Indirect  BEM 

The  essence  of  the  indirect  BEM  implementation  is  to  discretize  the  boundary  between 
cavity  fluid  and  surrounding  formation  into  a  set  of  small-sized  surfaces  called  elements. 
Each  element  is  a  ring-shaped  surface  with  height  dz  and  cavity  radius  r o-  Density  of  the 
fictitious  source  is  assumed  to  be  constant  on  each  element.  From  our  previous  boundary 
integral  formulation,  it  is  seen  that  a  fictitious  volume  source  for  the  fluid  and  a  fictitious 
source  vector  for  the  elastic  medium  are  needed  in  order  to  uniquely  describe  the  cavity 
source  radiation.  The  fictitious  source  vector  in  an  axially  symmetrical  system  consists 
only  of  the  vertical  and  radial  components.  Therefore,  we  have  on  each  element  i  three 
unknowns  to  determine  the  following:  a  fictitious  fluid  volume  source,  V/]  a  vertical  source 
for  elastic  medium,  and  a  radial  source  for  elastic  medium,  F[.  Our  goal  is  to  obtain 
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these  sources  on  each  element  so  that  displacement  fields  in  the  elastic  medium  or  the  fluid 


can  be  calculated. 

The  number  of  elements  is  restricted  by  the  computing  power  of  current  computers  and 
the  accuracy  specified.  In  practice,  three  elements  per  shortest  wavelength  are  used,  and 
the  number  of  elements  depends  on  the  time  window  and  the  fastest  wave  speed.  Thus,  the 


element  height  dz  satisfies 


min  (c, /?,  a) 


and  the  number  of  elements,  Ne,  is  given  by 

_  x  max(c,  /?,  a)  _  3/  x  tmax  x  ma x(c,/?,q) 
dz  min(c,  /3 ,  or) 

In  the  above,  the  frequency  and  maximum  time  to  be  calculated  for  are  denoted,  respec¬ 
tively,  by  /  and  tmax.  At  each  element,  the  usual  fluid/elastic  boundary  conditions  have 
to  be  satisfied:  i.e.,  continuity  of  normal  displacements  and  continuity  of  normal  stress  and 
vanishing  tangential  stress  of  the  solid, 


UrUi-UrU;  =  0 

arr\T=r+  ~  <7rrlr=r-  =  0 


0VZL_,+  =  0. 


These  boundary  conditions  are  satisfied  at  the  center  of  each  element.  Displacements  and 
stresses  (or  pressure)  at  the  center  of  each  element  result  from  contributions  from  fictitious 


sources  of  all  the  elements.  To  calculate  the  displacement  at  the  j-th  element,  due  to  the 
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source  at  the  i-th  element,  we  use  the  indirect  formulations  (1)  and  (3).  At  the  ji-th  element, 
the  boundary  conditions  become 

Ne  Nc  Nt 

=  D1 

1=1  i=i  i=i 

E  BjiVf  +  E  BJJ?  +  E  Brp;  =  D°"  (8) 

1  =  1  1=1  1=1 

E  q^+Eq,?;  =  D°". 

1=1  i=i 

In  the  above,  Al-{,  A"-  and  ATj{  represent  displacements  at  the  j-th  element  due  to  the  cylin¬ 
drical  volume,  and  the  vertical  and  radial  ring  sources  of  unit  strength  at  the  i-th  element, 
respectively.  They  are  integrals  of  the  Green’s  functions  with  respect  to  the  surfaces  of  cavity 
(bottom  +  wall).  The  B's  and  C’s  are  the  radial  and  tangential  stresses  at  the  j-th  element 
due  to  sources  at  the  i-th  element.  They  the  integrals  of  the  stress  Green’s  function.  The 
D's  are  exciting  fields  indicated  by  their  superscripts  at  the  j-th  element.  With  j  ranging 
from  1  to  Ne,  we  obtain  3  x  iVe  equations  that  can  be  easily  solved  for  the  3  x  Ne  unknowns. 
With  these  fictitious  source  densities  determined,  fields  inside  and  outside  of  the  cavity  can 
be  easily  obtained  using  equations  (1)  and  (3). 

Integration  of  the  Green’s  functions  is  done  in  two  steps.  First,  using  the  Sommerfeld 
integral  representation  for  function  e,kR/R,  we  transform  the  Green’s  functions  into  integrals 
with  respect  to  horizontal  wavenumber.  The  use  of  Sommerfeld  integral  representation  allows 
us  to  incorporate  vertical  layering  in  a  later  study.  The  r0dtp  part  of  surface  integration 
greatly  simplifies  the  integrand  of  the  wavenumber  integral  due  to  axial  symmetry.  Then  the 
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dz  part  (for  wall  elements)  and  the  dr0  part  (for  bottom  elements)  of  the  surface  integral  can 
be  analytically  evaluated.  The  integration  of  the  Green’s  function  is  reduced  to  a  horizontal 
wavenumber  integral.  Second,  this  wavenumber  integral  is  evaluated  using  the  discrete 
wavenumber  method  (Bouchon,  1981)  by  summation  over  the  horizontal  wavenumber. 

The  coefficient  matrix  of  the  equation  system  (8)  is  a  fully  populated  complex  matrix 
and  is  non-symmetrical.  This  is  often  regarded  as  the  disadvantage  of  BEM  as  compared  to 
FEM  or  FDM.  In  the  latter  two  methods,  tridiagonal  matrices  are  obtained  and  a  special 
faster  algorithm  exists  for  this  kind  of  matrix.  Nevertheless,  this  matrix  can  still  be  easily 
manipulated  as  the  number  of  elements  is  not  exceedingly  high  and  the  system  of  equations 
is  only  solved  once  for  each  frequency. 

Results 

The  situation  of  a  seismic  source  in  a  finite  cavity /borehole  arises  in  nuclear  explosions  in 
tunnels,  or  in  a  borehole  seismic  experiment  near  the  bottom  of  a  borehole.  Most  likely,  in 
the  first  case  the  tunnel  is  filled  by  air  (or  empty),  and  in  the  second  the  borehole  is  filled 
by  water.  In  either  case,  the  finiteness  of  the  borehole/cavity  will  play  a  role  in  the  far  field 
radiation  pattern  of  the  seismic  waves.  The  BEM  technique  outlined  above  can  be  used  to 
study  the  wave  radiation  from  this  special  source  geometry.  Some  steps  involved  in  the  BEM 
modeling  are  illustrated.  The  influence  of  a  finite  borehole  on  P  and  S  wave  radiation  is 
demonstrated. 
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For  a  finite  tunnel,  the  wall  of  the  tunnel  as  well  as  both  its  ends  have  to  be  discretized 
into  elements  (See  Figure  1).  On  each  of  these  elements,  we  need  to  determine  the  fictitious 
sources  on  the  solid  side  and  on  the  fluid  side  of  the  interface.  The  requirement  of  the 
three  boundary  conditions  leads  to  a  system  of  equations,  and  these  sources  can  thus  be 
determined.  In  a  homogeneous  medium,  this  BEM  technique  is  very  fast  due  to  the  following 
considerations: 

•  Let  Nu ,e  be  the  number  of  wall  elements.  Wall-to-wall  element  interactions  need  not 
be  calculated  for  each  of  the  N we  radiating  elements  and  each  of  the  Nwe  receiving 
elements,  which  would  involve  Nwe  x  Nwe  evaluation  of  a  wavenumber  integral.  Instead, 
calculating  the  interaction  of  one  receiving  element  with  all  Nwe  radiating  elements, 
which  incurs  only  Nwe  evaluation  of  integrals,  is  sufficient.  This  is  because  of  the 
permutation  and  symmetry  properties  of  the  displacement  and  stress  integrals  when 
radiating  and  receiving  elements  are  shifted  and  interchanged. 

•  If  we  call  one  end  of  the  tunnel  the  bottom,  and  the  other  the  top,  the  interactions 
between  each  top  element  and  each  bottom  element  and  their  self-interactions  can  be 
computed  as  efficiently  as  those  of  the  wall  elements. 

•  No  symmetry  property  can  be  used  for  computing  wall-to-top  element  interaction.  But, 
once  this  is  done,  the  wall-to-bottom  interaction  can  be  copied  from  the  wall-to-top 
calculations. 
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These  considerations  yield  a  tremendous  saving  of  computation  time. 

To  demonstrate  the  effect  of  the  finite  cavity  on  wave  radiation  from  an  explosion  source, 
two  situations  that  are  applicable  to  nuclear  and  non-nuclear  explosions  in  an  isotropic  ho¬ 
mogeneous  space  are  considered.  For  a  nuclear  explosion,  we  consider  that  the  compressional 
velocity  inside  the  cavity  is  the  shock  wave  velocity.  The  shock  wave  velocity  is  taken  as  10 
km/sec  (Lamb  et  ai,  1992).  For  a  non-nuclear  (i.e.,  chemical)  explosion,  the  compressional 
velocity  inside  the  cavity  is  330  m/sec  (sound  speed  in  air).  We  choose  an  isotropic  medium 
with  a  P  velocity  of  4.0  km/s,  a  S  velocity  of  2.2  km/s,  and  a  density  of  2.2  g/cm3  as 
a  homogeneous  background  model.  Four  tunnels  are  considered:  40  meters  long  with  a  5 
meter  radius,  40  m  long  with  a  10  m  radius,  80  m  long  with  a  5  m  radius,  and  80  m  long 
with  a  10  m  radius.  The  receivers  are  circularly  and  evenly  distributed  around  the  source 
at  a  distance  of  5  kilometers  from  the  center  of  the  tunnel  (see  Figure  1).  For  comparison 
with  the  radiation  patterns  of  sources  in  finite  cavity,  Figure  2  provides  the  P  and  S  wave 
radiation  patterns  for  an  explosion  source  in  an  air-filled  infinite  cylindrical  cavity. 

Figures  3  and  4  show  the  radial  and  tangential  components  for  the  explosion  source  placed 
at  1  and  18  meters,  respectively,  from  one  end  of  the  tunnel  with  a  5  meter  radius  and  40 
meter  length.  In  this  case,  the  compressional  velocity  inside  the  tunnel  is  chosen  as  air.  The 
center  frequency  is  1  Hz.  These  figures  display  the  seismograms  and  the  radiation  pattern 
of  the  waves  simultaneously.  The  radial  seismograms  show  the  radial  particle  motion  and 
hence  the  P  wave  radiation.  The  tangential  seismograms  show  the  transverse  motion,  i.e., 
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the  shear  wave  radiation.  Figure  5  shows  the  seismograms  from  a  nuclear  explosion  (shock 
wave  velocity  inside  cavity)  located  18  meters  from  one  end  of  the  tunnel.  The  maximum 
displacement  of  the  P  and  5  wave  at  the  receivers  for  two  center  frequencies  (1  and  2  Hz)  is 
plotted  in  Figures  6  through  13  for  nuclear  and  non-nuclear  explosions  at  various  positions 
along  the  tunnel  axis. 

Figures  3  through  13  show  that,  when  an  explosion  source  is  at  the  center  of  a  finite 
cylindrical  tunnel,  its  S  wave  radiation  pattern  is  basically  the  same  as  that  of  an  explosion 
source  in  an  infinite  borehole  (Heelan,  1953;  Lee  and  Balch,  1982;  Meredith,  1990).  The 
P  wave  radiation  pattern,  on  the  other  hand,  is  close  to  that  of  a  force  along  the  cavity 
axis.  Slightly  away  from  the  center  of  the  cylindrical  cavity,  the  S  radiation  pattern  of  a 
non-nuclear  explosion  (air  velocity)  changes  drastically.  Near  the  end  of  the  tunnel,  S  and  P 
wave  radiation  patterns  from  a  non-nuclear  explosion  are  essentially  those  of  a  single  force 
along  the  tunnel  axis,  and  the  radiated  waves  have  their  maximum  amplitude  when  the 
source  is  close  to  the  tunnel  end.  In  contrast,  the  nuclear  explosion  (shock  wave  velocity 
inside  the  cavity)  produces  similar  radiation  patterns  for  different  source  positions  along 
the  axis  of  the  cylinder.  Shear  waves  are  generated  at  each  position,  with  the  shear  wave 
amplitude  varying  slightly  with  the  position  of  the  source  especially  at  the  end  of  the  cavity. 
The  dipole  source  pattern  of  shear  wave  radiation  is  very  consistent  for  all  the  examples. 

Both  ends  of  the  tunnel  are  responsible  for  these  peculiar  radiation  patterns.  Energy 
from  an  explosion  is  mostly  trapped  as  tube  waves.  The  impact  of  the  direct  wave  and 
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the  tube  wave  on  both  ends  of  the  cavity  is  much  stronger  than  on  the  cavity  wail.  When 
the  explosion  source  is  at  the  center,  both  ends  experience  the  same  pressure  at  the  same 
time,  and  in  the  far  field  they  act  like  a  single  force  couple.  This  results  in  an  S  wave 
radiation  pattern  similar  to  that  of  an  infinite  borehole.  The  P  wave  pattern  results  from 
the  fact  that  receivers  close  to  one  end  of  the  cavity  detect  more  P  wave  energy  from  this 
end  than  from  the  other.  As  the  source  moves  away  from  the  center,  one  end  of  the  tunnel 
dominates,  resulting  in  the  kinds  of  radiation  patterns  observed  especially  for  a  non-nuclear 
explosion.  For  nuclear  explosions,  the  high-speed  propagation  velocity  inside  the  cavity  could 
not  produce  the  asymmetric  results  seen  in  the  case  of  non-nuclear  sources.  In  this  case,  the 
entire  tunnel  surface  radiates  energy  at  the  same  time  as  a  volume  source. 

Conclusions 

A  general  numerical  algorithm  based  on  the  boundary  element  method  has  been  used  to 
model  seismic  waves  radiated  from  explosion  sources  inside  finite,  cylindrical  tunnels.  This 
method  has  the  capability  to  calculate  both  near-field  and  far-field  radiation  and  the  scale 
difference  is  properly  dealt  with.  We  computed  source  radiation  patterns  for  two  different 
explosion  source  types  placed  along  the  axis  of  a  finite  tunnel.  The  radiation  patterns  from 
non-nuclear  sources  depend  strongly  on  the  position  of  the  source  along  the  axis.  Sources 
near  the  tunnel  wall  produce  the  highest  amplitudes  and  behave  like  a  simple  point  force. 
Radiation  patterns  from  nuclear  sources  do  not  vary  significantly  with  source  position  owing 
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to  the  very  high  compressional  velocity  inside  the  tunnel.  All  the  sources  and  source  positions 

considered  produced  significant  shear  waves,  showing  that  a  finite  cylindrical  cavity  provides 

a  mechanism  for  shear  wave  generation  from  explosions. 
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FINITE  TUNNEL 


a  =  4  km/sec 
(3  =  2.2  km/sec 


Figure  1:  Schematic  diagram  of  the  element  of  the  tunnel  and  the  source  and  receivers 
configurations. 


18 


Infinite  Cylindrical  Cavity 


Figure  2:  Analytical  P  ( solid  line )  and  SV  ( dashed  line )  wave  radiation  pattern  for  an 
explosion  source  at  the  center  in  an  infinite  air-filled  cylindrical  cavity. 
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Radial 


Tangential 


Radius  =  5.  meter 
Length  =  40.  meter 
Zs  =  1.  meter 
Rr  =  5000.  meter 
f_c  —  1.  Hz 

A_max  =  0.1 52077 13E-09 


Figure  3:  Radial  and  tangential  seismograms  at  a  distance  of  5  km  from  a  finite  tunnel  with 
a  length  of  40  meters  and  a  5  meter  radius.  The  center  frequency  is  1  Hz.  The  source  is 
placed  1  meter  from  one  end  of  the  tunnel.  The  compressional  velocity  inside  the  tunnel 
is  330  m/sec  (air). 


Figure  4:  Same  as  described  in  Figure  3  but  the  source  is  18  meters  from  one  end  of  the 


Radius  =  5  Radius  =  10 


Cavity  length  =  40  meter  (air,  1Hz) 


Figure  6:  Radiation  pattern  of  P  and  S  waves  from  a  finite  tunnel  of  length  40  meters  for 
different  source  positions  along  the  tunnel  with  compressional  velocity  330  m/sec  ( non¬ 
nuclear  source).  The  center  frequency  of  the  waveforms  is  1  Hz.  The  left  column  shows 
the  results  for  the  tunnel  with  a  5  meter  radius  and  the  right  column  is  for  a  10  meter 
radius.  The  source  position  from  one  end  of  the  tunnel  is  also  shown  in  the  middle 
column.  The  amplitudes  are  plotted  on  the  same  scale,  but  the  multiplication  factor 
(e.g.,  x3  for  radius  5  and  source  position  1  meter  from  the  end  of  the  tunnel)  is  also 
shown  where  required. 


Radius  =  5 


Radius  =  10 


Cavity  length  =  40  meter  (Shock  wave,  1  Hz) 


Figure  10:  Same  as  described  in  Figure  6.  The  compressional  velocity  inside  the  tunnel 

10  km/sec  as  shock  wave  (nuclear  explosion). 
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Radius  -  5 


Radius  =  10 


Figure  12:  Same  as  described  in  Figure  10.  The  center  frequency  is  2  Hz. 
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Cavity  length  =  80  meter  (Shock  wave,  4  Hz) 


Same  as  described  in  Figure  11.  The  center  frequency  is  2  Hz. 


EXPERIMENTAL  STUDY  OF  SCATTERING  FROM  A 


HIGHLY  IRREGULAR,  ACOUSTIC-ELASTIC  INTERFACE 
Summary 

In  this  study,  we  experimentally  and  numerically  investigate  the  scattering  of  an  acoustic  P 
wave  incident  on  a  highly  irregular,  random  acoustic-elastic  interface  to  determine  whether 
enhanced  backscattering,  already  identified  numerically  for  SH  and  P-SV  waves,  occurs.  Nu¬ 
merically,  the  problem  is  solved  by  coupling  the  Somigliana  identity  for  an  elastic  medium 
with  Green’s  second  integral  theorem  for  pressure  in  a  fluid.  Exact  integral  expressions  for 
the  scattered  pressure  in  the  acoustic  medium  are  then  obtained.  Experimentally,  a  glass 
etching  process  using  photoresist  templates  with  Gaussian  statistics  allowed  for  the  genera¬ 
tion  of  a  characterized  random  interface.  This  3-D  interface  has  approximately  a  Gaussian 
correlation  function  and  a  Gaussian  height  distribution.  A  method  was  also  developed  by 
which  identical  interface  geometries  with  differing  material  contrasts  can  be  physically  cre¬ 
ated.  This  approach  involved  using  the  glass  surface  as  an  epoxy  mold.  Experiments  were 
carried  out  on  the  glass  surface  in  M.I.T.’s  Earth  Resources  Laboratory’s  ultrasonic  labora¬ 
tory.  Two-dimensional  numerical  results  predict  the  3-D  experimental  results  well  at  small 
incident  angles.  Both  numerical  and  experimental  results  strongly  support  the  presence  of 
"enhanced  backscattering.”  However,  more  experiments  are  required  to  fully  constrain  the 
properties  of  the  retroreflective  peak.  Finally,  fundamental  differences  between  2-D  and  3-D 
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scattering  mechanisms  appear  to  form  at  larger  incident  angles. 


Introduction 

In  laboratory  experimentation  a  lack  of  control  over  the  statistical  parameters  of  a  given 
random  model  can  easily  make  experimental  results  ambiguous.  In  the  case  of  irregular 
interfaces,  the  height  probability  distribution  and  the  correlation  lengths  of  the  interface 
may  not  be  well  constrained,  the  interface  may  not  be  stationary  in  space,  and  the  interface 
may  contain  a  wide  variety  of  length  scales.  Each  of  these  experimental  uncertainties  makes 
comparisons  with  numerical  models  difficult,  if  not  impossible.  It  is  the  goal  of  this  project 
to  physically  fabricate  a  random  interface  which  is  stationary  in  space  with  both  a  simple 
height  probability  distribution  and  a  simple  transverse  correlation  function  so  that  these 
experimental  results  can  be  easily  compared  with  the  corresponding  numerical  results. 

Experimentally,  the  accurate  generation  of  a  Gaussian  surface  is  very  important  since 
Gaussian  interfaces  are  mathematically  convenient  and  widely  used  in  scattering  studies. 
Many  theoretical  formulations  in  the  literature  apply  the  simple  properties  of  a  Gaussian 
correlation  function  to  random  surfaces.  Examples  can  be  found  in  Prange  and  Toksoz 
(1990),  KnopofF  and  Hudson  (1964,  1967),  Haddon  (1978),  and  Kuperman  and  Schmidt 
(1989).  Exponential  correlation  functions  have  also  been  used  extensively  (e.g,  Wu  and  Aki, 
1985;  Frankel  and  Clayton,  1986)  and  in  many  instances  give  a  good  description  of  field 
observations. 
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In  this  study,  we  choose  to  fabricate  and  model  an  interface  with  a  Gaussian  correla¬ 
tion  function.  The  statistical  parameters  of  the  interface  are  chosen  so  that  the  incident 
wavelength  has  the  same  length  scale  as  the  correlation  length  of  the  irregularities.  In  ad¬ 
dition,  the  average  slope  of  the  interface  is  large  and  the  approximate  techniques  such  as 
the  Kirchoff,  Born,  and  the  geometrical  ray  approach,  break  down  into  multiple  scattering 
mechanisms  such  as  “enhanced  backscattering”  and  “shadowing,”  and  play  strong  roles  in 
wave  scattering. 

What  is  enhanced  backscattering  from  an  acoustic-elastic  halfspace?  By  definition,  en¬ 
hanced  backscattering  or  ‘retroreflectance’  is  the  enhancement  of  energy  scattered  back  in  the 
direction  of  the  source.  O’Donnell  and  Mendez  (1987)  were  the  first  to  propose  the  hypoth¬ 
esis  that  time-reversed  paths  were  responsible  for  enhanced  backscattering.  This  hypothesis 
was  further  strengthened  by  Maradudin  et  al.  (1991)  who  showed  that  retroreflectance  for 
energy  double  scattered  from  the  interface.  More  support  came  from  Schultz  and  Toksoz 
(1993b),  who  showed  that  full  elastic  seismic  scattering  is  consistent  with  this  hypothesis 
since  enhancement  is  observed  on  P-to-P  and  S-to-S  scattering  and  not  on  P-to-S  and  S-to-P 
scattering.  The  extension  of  time-reversed  paths  is  easily  extended  to  the  acoustic-elastic 
case.  Take  for  instance  the  peak- valley  sequence  shown  in  Figure  1.  If  an  incident  P-wave, 
shown  by  the  solid  line,  diffracts  from  point  1,  it  will  propagate  as  a  P-wave  to  point  2 
and  then  diffract  at  some  angle  into  the  upper  medium  again  as  a  ±  -  wave.  For  most  waves 
travelling  away  from  the  interface,  enhancement  will  not  occur.  However,  if  the  diffracted 
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wave  travels  directly  back  towards  the  source,  an  incident  P-wave  can  be  found  travelling 
exactly  in  the  reverse  path — propagating  from  point  2  to  point  1  and  again  travelling  back 
towards  the  source  as  shown  by  the  dashed  line.  In  this  case,  the  time-reversed  path  in¬ 
terferes  constructively  with  the  forward  path  and  contributes  additional  energy  towards  the 
source,  resulting  in  enhanced  backscattering.  Using  a  simple  phase  argument,  some  proper¬ 
ties  of  enhanced  backscattering  can  be  derived.  The  peak  width  can  be  written  as  A 6,  = 
where  Ad,  is  the  angular  width  of  the  peak,  A  is  the  incident  wavelength,  and  /  is  the  mean 
free  path  of  the  interface  or,  in  other  words,  the  average  distance  a  wave  propagates  between 
points  1  and  2  along  the  interface. 

Using  this  phase  approach  other  path  geometries  may  also  contribute  to  enhanced  backscat¬ 
tering.  For  example,  if  a  wave  encounters  the  interface  and  multiply  scatters  three  times  as  a 
P-wave  to  send  energy  back  in  the  direction  of  the  source,  then  a  time-reversed  path  may  be 
found  that  also  sends  energy  back  towards  the  source.  In  the  same  manner,  many  multiply 
scattered  paths  sending  energy  back  towards  the  source  can  be  found  in  the  acoustic-elastic 
case.  However,  as  a  result  of  energy  loss  with  each  diffraction  from  the  interface  due  to  both 
transmission  through  the  interface  and  additional  spreading,  it  seems  reasonable  that  the 
double-scattered  path  contributes  the  majority  of  retroreflective  energy. 

In  this  paper,  we  first  briefly  summarize  the  numerical  formulation  used  to  model  scatter¬ 
ing  from  a  randomly  irregular  acoustic-elastic  interface.  This  Somigliana  identity  approach 
is  based  on  the  work  of  Schultz  and  Toksoz  (1993a, b).  Next,  the  construction  of  the  ran- 
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dom  glass  surface  is  discussed  and  the  ultrasonic  apparatus  for  measuring  the  scattering 
properties  of  the  interface  is  then  described.  The  third  section  compares  the  experimental 
scattering  results  obtained  in  our  in-house  ultrasonic  water  tank  with  the  numerical  mod¬ 
els.  The  scattering  properties  are  discussed  in  detail  along  with  a  discussion  of  probable 
scattering  mechanisms.  Finally,  we  discuss  the  differences  between  the  2-D  synthetic  results 
and  the  3-D  experimental  data.  We  then  propose  possible  differences  between  2-D  and  3-D 
scattering  mechanisms. 

Theory 

The  numerical  approach  described  in  this  section  and  the  notation  used,  follows  that  of 
Schultz  and  Toksoz  ( 1993a, b).  Since  the  approach  here  is  very  similar  to  the  SH  and  P-SV 
cases,  we  give  only  a  brief  outline  of  the  theoretical  approach.  We  first  express  the  total 
scattered  displacement  at  any  point  within  two  volumes  of  elastic  material  exactly  with  the 
Somigliana  representation  theorem  (e.g..  Aki  and  Richards,  1980).  Simplifying  this  theorem 
to  a  2-D  case  gives  a  set  of  four  integral  equations 

=  Jv4l>(n)G^-,v)dV(1) 

-(—I)'  /  <tf(x'){(cW,fti(x')aG"(x ;x')/3z>!'>(x')  (9) 

J  s 
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where  gradients  are  zero  in  the  x2-direction.  Following  the  notation  of  Schultz  and  Toksoz 
(1993b),  T^(x)  is  the  traction  vector  along  the  interface  in  the  fluid  (/  =  f)  and  the  solid 
(/  =  s),  and  we  have  assumed  all  surfaces  are  far  enough  away  that  only  the  surface,  S(x), 
separating  the  two  volumes,  contributes  to  the  final  displacement.  Gnp(x;x  )  is  a  Green’s 
function  giving  the  n-component  of  displacement  at  x  resulting  from  a  point  force  in  the 
p-direction  at  x  ,  clJpq  is  the  elasticity  tensor,  and  H[ i]  is  a  function  that  takes  a  value  of 
0  or  1  depending  on  whether  the  point  x  lies  outside  or  inside  the  volume  of  interest,  i, 
respectively.  We  assume  that  the  media  are  homogeneous  and  isotropic,  so  the  constitutive 
relation  can  be  written  with  the  aid  of  the  elasticity  tensor  as 

ru(x)  =  AO(x)6;j  +  p(ui,j(x)  +  uiii(x)],  (10) 

where  0(x)  =  Uk,k(%)  is  the  dilatational  parameter  and  A  and  /z  are  Lame  parameters  for 
the  medium. 

In  this  work  the  upper  medium  is  acoustic,  supporting  propagation  of  only  dilatational 
waves  while  the  lower  medium  is  taken  as  elastic.  The  boundary  separating  these  two  media 
is  shown  in  Figure  2.  The  boundary  conditions  for  the  resulting  acoustic-elastic  boundary, 
can  be  written  in  the  general  form 

0  •  yl/)(x)r3=c(ri)  =  0  •  u(3,(x)X3=c(ll), 

T(/)(x)x,=c(x,)  =  T(i)(x)x3=c(x1),  (11) 

!1  x  I  (s)(xL3 =<(*,)  =  0, 
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where  the  surface  profile  function  is  taken  to  be  13  =  C(xi)  and  t^e  unit  normal  vector  along 
the  interface  can  be  expressed  as 

n  =  (-C/(x1))l)[l  +  (C'(a;1))2]^.  (12) 

The  first  boundary  condition  represents  the  continuity  of  normal  displacement  and  the  latter 
two  conditions  together  represent  the  continuity  of  normal  stress.  Referring  to  the  require¬ 
ment  for  continuity  of  normal  displacement  one  can  expand  the  first  term  of  the  volume 
integral  (9)  as 

n,c|/p)gti!/)(x)  =  XU)ntu\n(x)6pq,  (13) 

which  follows  from  the  lack  of  rigidity  in  the  acoustic  medium,  =  0.  The  last  two 
boundary  conditions  in  (11)  infer  the  continuity  of  fluid  pressure  at  the  interface.  Upon 
combining  these  two  boundary  conditions  with  the  constitutive  relation,  (10),  the  traction 
in  the  solid  can  be  expressed  at  the  interface  as 

T^s)(x)  =  SU){x)np  =  A  (14) 

where  5^*(x)  is  the  fluid  pressure  in  the  fluid.  Finally,  referring  to  eq.  (13)  and  comparing  it 
with  a  similar  expansion  in  the  elastic  medium  it  is  clear  that  the  equality  u-^(x)  =  £/,(x)  = 
u-^(x)  implies  that  the  normal  displacement  is  continuous,  or  nku[^(x)  =  nfcuj^(x)  = 

nki'k(x)- 

Taking  our  volume  of  interest  to  be  the  upper  acoustic  medium,  placing  the  incident  wave 
in  the  acoustic  medium,  and  substituting  the  final  form  of  the  boundary  conditions  (11),  the 
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integral  equation  in  the  lower  elastic  medium  can  be  written  as 


(15) 


where  the  unknowns  are  the  fluid  pressure,  S’^(x),  and  the  displacement,  U(x),  along  the 
interface.  S^(x)  should  not  be  confused  with  the  surface  function  referred  to  previously. 

In  the  acoustic  medium,  the  surface  integral  can  be  greatly  simplified.  Since  shear  waves 
can  not  propagate  in  the  acoustic  medium,  one  of  the  two  integral  equations  (9)  in  the 
upper  medium  is  redundant  and  can  be  combined  into  one  equation.  We  note  that  the  fluid 
pressure  in  the  acoustic  medium,  S^\x),  satisfies  the  basic  wave  equation 


V25^(x)  +  (1 k/)25(/)(x)  =  0, 


(16) 


assuming  no  sources  are  present  in  the  medium.  The  Helmholtz  potential  for  the  displace¬ 
ment,  <^U)(> c),  also  satisfies  a  similar  wave  equation.  Transforming  to  the  frequency  domain 
and  differentiating  we  find  that  d^*(x)  is  directly  related  to  the  P  wave  pressure, 

o(/)(x)  =  -(pifW)-lSu){x),  (17) 


Using  this  relation,  the  normal  displacement  at  the  interface  can  be  expressed  as 

— =  -P{f)u2[nku[*\x))  =  —A  k{Tf)2{nku[f\x))  (18) 

and  using  the  Green’s  function  for  the  pressure,  G^\x-,x'),  the  integral  equation  in  the 
acoustic  medium  can  be  expressed  with  the  aid  of  Green’s  second  integral  theorem  as 

5(/>(x)  =  .s'(/)(x)mcld  (19) 
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+  [+°°  dx[[\WSu\x')dGU +  A</>a4/,2G(/)(x;  s,)nfctil/)(jt,)]> 

J —oo  UTL 

where  we  have  utilized  (18).  Equations  (15)  and  (19)  together  consist  of  three  integral 
equations  with  three  unknown  functions.  These  integral  equations  now  express  the  total 
scattered  field  in  both  media.  Letting  x3  — *  ^+(x1),  the  final  set  of  coupled  integral  equations 


can  be  written  as, 


S(x)  =  S^(x)ineid 


/-f-OO 

dx\lS(i‘)Tu)(x\s)  - 

-oo 


in  the  acoustic  medium,  and 


0  =  -  r“,<fx;[f/,(x')rl‘l(5ls')-^'-|(x|x')S(x')].  (21) 

J  -OO  f 

in  the  elastic  medium.  VVe  have  defined 

r«'l(x|x')  =  V/,aG,2^')U-c,„.  (22) 

D^xlx')  =  -A,'ity'r-6',,)(x;x')„,|„=ailJ,  (23) 

rM(x\x')  =  7r(‘l(si5')i.J=<(.,), 

^’(xlx)  = 

in  the  respective  media.  Now  the  unknown  source  strength  functions,  which  we  eventually 
solve  for.  can  be  expressed  as  a  function  of  x\  alone 

S(xr)  =  0{/)(x)U=«x,),  (24) 


L^i(xi)  —  Li  (x)|x3=c(xi)) 
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where  we  have  normalized  the  pressure  term  with  respect  to  to  ensure  that  the  final 
numerical  conditions  are  well  conditioned. 

The  scattered  field  in  the  acoustic  medium  can  now  be  expressed  completely  in  terms  of 
the  unknown  source  functions,  (24).  The  cartesian  coordinate  form  of  the  Green’s  function 
in  the  fluid  can  be  written  as 


G(/)(x;x')  = 


f+OO  p>'i(*l-*i)+«'*4/>lX3-:r3l 


47rA(^) 


where 


ldf)  =  ((k{Tf))2-k2)>,  Im(k(3f))  >  0, 

which  corresponds  to  a  pressure  source,  P(x)  =  —  —  x'3)  applied  at  x* 

in  the  fluid.  Substituting  (25)  into  the  surface  integral  (19)  the  scattered  field  at  any  point 
x3  >  C(xi)maI  in  the  fluid  can  be  decomposed  into  a  summation  of  plane  waves 


S{x)  scat  =  +  J 


+°°  dk 

•°°  2  vk&n 


Rf(kuj)e 


+t/cxi  +ikjX3 


where  the  amplitude  coefficient  takes  the  form 


ff/(M  =  5  I*™  <lx[{iS(x')(kC(x')  -  k'") 

1  J —oo 


-  D3(x'))]e 


To  reduce  the  computational  demand  of  this  approach,  the  incident  wave  is  expressed  as 
a  narrow  Gaussian  beam  source  following  Maradudin  (1990).  This  allows  for  a  reduction 
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in  the  length  of  integration  along  the  acoustic  boundary,  since  only  a  small  portion  of  the 
interface  is  excited  by  the  incident  beam.  The  pressure  of  a  Gaussian  beam  incident  at  an 
angle,  0O>  in  the  acoustic  medium  can  be  written  as 


5(/)(x), 


_  i*nfio-r3CO*®o)[l  +  M'(X)]e-((iicose0+*3»*"®o)/«')2 


(28) 


where 

IV(x)  =  — {-^(xxcosQo  +  x3sin90)2  -  1],  (29) 

k[Tniu2  w2 

which  is  an  approximation  to  the  wave  equation,  (16),  and  is  valid  so  long  as  ^  >>  1, 
where  w  is  the  half-width  of  the  Gaussian  beam.  We  also  express  w  =  hcos00,  where  h  is 
the  half-width  of  the  incident  beam  as  seen  on  the  plane  x3  =  0. 

An  approximation  to  the  Fourier  reflection  coefficient  can  now  be  written  in  terms  of 
the  amplitude  coefficient,  (27).  Normalizing  this  amplitude  coefficient  by  the  amplitude  of 
the  incident  plane  wave  having  a  wave  vector  corresponding  to  the  angle  of  incidence,  the 
reflection  coefficient  can  be  expressed  as 


where 


R{ku)  = 


MM1  , 

'ly/ir  kj\u 


(30) 


?(^s)  =  /  (lx[[i S (x[) kj  \sm 0 ,(( x[)  ~  cos  9 s)  (31) 

■/+» 

-4/)2(C'(*i)0i(*i)  -  ^3(®i))]c"i^)(*infl**,,+c“tf*<'(r,,))» 


liich  is  comparable  in  amplitude  to  the  Fourier  reflection  coefficient  calculated  for  a  single 
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incident  plane  wave,  and  we  let  k  =  k^  sin  0,  and  k  =  cos6a.  Note  that  this  normaliza¬ 
tion  is  different  from  the  Differential  Reflection  Coefficient  of  Schultz  and  Toksoz  ( 1993a, b). 

Appendix  A  describes  how  these  integral  equations  can  be  solved  numerically  following 
the  approach  of  Schultz  and  Toksoz  ( 1993a, b).  We  show  that  the  solution  to  this  acoustic- 
elastic  case  can  be  expressed  completely  as  a  combination  of  the  Green’s  functions  for  the 
P-SV  and  SH  case,  where  the  shear  velocity  of  the  SH  wave  is  changed  to  the  P-wave 
velocity  of  the  acoustic  medium,  so  as  to  reflect  the  acoustic  Green’s  function.  The  final 
coupled  integral  equations  are  then  transformed  to  a  coupled  set  of  matrix  equations  and 
solved  using  LU  decomposition. 

Numerical  Analysis 

In  this  paper  two  types  of  interfaces  are  modeled  numerically.  Both  interfaces  have  a  Gaus¬ 
sian  distribution  about  the  mean,  where  S2  =  (£2(xi))  is  the  mean-square  departure  of  the 
surface  from  flatness.  The  first  interface  studied  has  a  correlation  function 

w(|x,-x;i)  =  r!(c(x.K(*;»,  (32) 

described  by  a  Gaussian  correlation  function,  W(|xi|)  =  exp(— x\/a2).  The  correlation 
length  for  a  Gaussian  interface  is  approximately  equal  to  the  average  distance  between  ad¬ 
jacent  peaks  and  valleys  along  the  interface.  The  interface  can  also  be  described  in  terms 
of  the  rms  slope  of  surface,  <p,  which  we  will  refer  to  often.  This  rms  slope  can  be  written 
as  o  =  tan~l(^L).  The  second  surface  studied  has  an  exponential  correlation  function, 
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lV(|xi|)  =  exp(—xi/a). 

Averaging  over  an  ensemble  of  realizations  of  these  surfaces,  we  display  the  final  scattered 
pressure  as  a  statistical  average  that  follows  the  approach  of  Schultz  and  Toksoz  ( 1993a, b). 
The  statistical  characteristics  of  the  scattered  pressure  can  then  be  analyzed  for  possible 
ttering  mechanisms.  Although  the  coherent  and  incoherent  properties  of  the  scattered 
>d  were  computed  for  seismic  analysis,  only  the  total  mean  squared  contribution  to  the 
Reflection  Coefficient  (RC)  is  presented  in  the  following  sections.  The  total  mean  squared 
contribution  to  the  RC  can  be  written  as 

=  wh{Km2)-  (33) 

This  gives  the  average  squared  pressure  reflected  into  the  upper  medium  as  a  function  of  the 
scattering  angle,  6S,  given  one  incident  beam  angle,  0O.  The  square  root  of  this  RC  is  used 
for  comparison  with  experimentally  recorded  amplitudes. 

Experimental  Procedure 

The  experimental  procedure  involved  submerging  a  solid  elastic  model,  in  this  case  a  glass 
block,  into  an  experimental  water  tank,  essentially  creating  an  acoustic-elastic  interface 
at  the  boundary  between  the  block  and  the  surrounding  water.  The  first  portion  of  the 
experiment  entailed  generating  a  characterized  random  surface  that,  if  successful,  would 
have  predetermined  Gaussian  statistics.  The  second  portion  involved  constructing  a  motor 
driven  measurement  device  that  could  accurately  (to  within  a  fraction  of  a  degree)  measure 
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various  realizations  of  the  fluid-glass  boundary. 


Random  Interface  Generation 

Fabricating  the  randomly  irregular  scattering  surface  was  the  most  challenging  part  of  this 
project.  Numerous  irregular  surfaces  were  designed.  Models  ranged  from  irregular  distri¬ 
butions  of  glass  beads  to  roughened  cement  surfaces.  In  addition,  random  distributions  of 
gravel  held  together  by  epoxy  were  tested  along  with  naturally  rough  granite  and  sandstone 
surfaces.  Unfortunately,  these  models  either  did  not  give  proper  control  over  statistical  pa¬ 
rameters  or  were  extremely  heterogeneous  at  the  ultrasonic  level,  making  comparisons  with 
numerical  models  very  difficult.  After  much  experimentation  the  most  promising  approach 
became  the  fabrication  of  a  random  glass  surface  using  an  etching  procedure. 

The  irregular  glass  surface  was  designed  using  a  solid  glass  block  and  a  standard  etching 
process.  First  the  cylindrical  glass  block  shown  in  Figure  3  was  cast  using  a  graphite  mold. 
After  one  week  of  annealing  the  block  had  a  final  diameter  of  19.5  cm  and  a  height  of  7.5  cm. 
The  upper  surface,  that  was  exposed  to  air,  contracted  slightly  due  to  the  high  expansion 
coefficient  of  Na  glass  and  resulted  in  a  slightly  concave  surface.  Therefore,  the  lower  plane 
surface  of  the  block  was  etched. 

The  general  theory  behind  etching  a  specific  surface  is  shown  in  Figure  4.  Take,  for 
instance,  the  fabrication  of  the  valley  shown  in  the  upper  left  portion  of  the  figure.  In  this 
case,  the  valley  is  divided  into  a  number  of  discrete  depth  intervals.  Photoresist  templates 
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are  designed  to  match  the  geometry  of  the  valley  at  each  discrete  depth.  The  first  template 
is  glued  to  the  smooth  glass  surface  that,  in  the  case  of  a  valley,  contains  only  a  small  gap. 
The  surface  is  then  exposed  to  a  constant  pressure  of  high  velocity  particles  which  chip  any 
portion  of  the  surface  not  covered  by  photoresist.  After  a  set  time  which  depends  on  the 
compressed  air  pressure  and  the  sand/air  mixture  ratio,  the  exposed  glass  is  etched  to  a 
target  depth.  The  first  layer  of  photoresist  is  then  removed,  and  the  next  layer  is  attached 
exposing  a  larger  portion  of  the  surface  to  the  incoming  sand  particles.  The  glass  surface  is 
exposed  again  to  sand  particles  for  the  same  amount  of  time.  The  valley  is  now  broader  and 
twice  as  deep.  Adding  each  template  in  a  similar  manner,  the  desired  valley  is  etched  into 
the  glass  surface. 

To  achieve  the  desired  random  interface,  the  Gaussian  surface  described  in  the  previous 
section  was  first  numerically  generated.  Both  the  transverse  correlation  length,  a,  and  the 
standard  deviation  of  the  height,  6,  of  the  interface  were  specified  as  1  mm  and  .71  mm, 
respectively,  giving  an  rms  slope  of  45°.  After  generation,  the  Gaussian  surface  was  dis¬ 
cretized  into  six  individual  depth  levels,  with  each  level’s  thickness  equal  to  one  standard 
deviation  of  the  surface.  The  templates  shown  in  Figure  5  were  successively  glued  to  the 
surface  and  each  template  was  exposed  to  high  velocity  sand  particles  normally  incident  on 
the  surface  for  approximately  360  s.  A  simple  lateral  sweeping  motion  of  the  sand  blaster 
was  used  to  cover  the  whole  template.  The  blaster  operated  at  a  pressure  of  125  kPa  (as 
IS  psi)  with  the  glass  surface  0.3  m  from  the  nozzle  of  the  sand  blaster.  In  general,  the 
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correlation  length  of  the  surface  was  controlled  by  the  template  design,  and  the  standard 
deviation  of  the  interface  was  controlled  by  the  time  that  each  template  was  exposed  to  the 
sand  blast. 

As  the  etching  proceeded,  we  found  that  the  standard  deviation  and  correlation  length 
were  difficult  to  control  precisely.  We  observed  during  the  sand  blasting  process  that  particles 
impacting  the  surface  tended  to  scatter  many  times  within  valleys,  generally  broadening  the 
valley  width.  In  addition,  small  narrow  peaks  tended  to  chip  away  far  faster  once  completely 
exposed  to  the  sand  blast,  removing  the  linearity  assumed  in  the  design  of  the  surface.  As 
we  show  below,  this  results  in  a  longer  correlation  length  then  expected.  In  this  study,  the 
target  ‘rms’  slope  of  the  interface  was  45°  and  the  desired  correlation  length  was  1  mm. 

The  Scattering  Instrumentation 

Once  the  irregular  glass  surface  was  created,  an  automated  scattering  apparatus  was  used 
to  measure  the  scattering  properties  of  the  interface.  The  two  different  flat-bottomed  trans¬ 
ducers  used  to  create  a  beam  source  were  a  Panametrics  12.7  mm  diameter  transducer  (1.5 
MHz,  A  =  1.0  mm  in  water)  and  a  Panametrics  25.4  mm  diameter  transducer  (0.5  MHz, 
A  =  3.0  mm  in  water).  The  detectors,  which  were  also  Panametrics  flat-bottomed  trans¬ 
ducers,  consisted  of  a  6.4  mm  diameter  transducer  (1.5  MHz,  A  =  1.0  mm  in  water)  and  a 
12.7  mm  diameter  transducer  (0.5  MHz,  A  =  3.0  mm  in  water),  respectively.  The  detectors 
were  chosen  such  that  they  were  sensitive  only  to  waves  approaching  nearly  perpendicular 
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to  the  surface  of  the  transducer,  limiting  the  energy  recorded  to  waves  which  approach  in 
line  with  the  transducer  axis.  Given  the  source  parameters,  the  resulting  source  radiation 
pattern  was  a  beam  of  energy  with  the  majority  of  the  source  energy  travelling  in  the  forward 
direction.  As  shown  in  Figure  6  the  source  radiation  patterns  show  some  spreading  of  the 
beam  although  further  tests  showed  that  this  slight  spreading  did  not  significantly  affect  the 
results. 

Experimental  data  was  recorded  in  our  in-house  water  tank,  described  in  Appendix  B. 
The  experimental  geometry  used  to  measure  the  surface  scattering  is  shown  in  Figure  7. 
The  glass  block  was  located  at  the  center  of  the  experiment  and  the  source  and  detector  were 
then  stepped  in  a  semicircle  about  an  axis  of  rotation  running  parallel  to  the  general  trend 
of  the  irregular  fluid-glass  interface.  In  each  experiment  the  source  was  placed  at  a  constant 
incident  angle,  Oq,  and  a  constant  .35  m  distance  from  the  axis  of  rotation.  The  recording 
angle  was  then  controlled  by  mounting  the  detector  on  a  motor-driven,  rotating  arm  that 
held  the  detector  .30  m  from  the  given  axis  of  rotation.  The  arm  was  then  rotated  in  0.9° 
steps  about  this  axis  of  rotation.  Therefore,  the  recorded  energy  represents  scattering  in  the 
plane  of  incidence.  As  a  result,  it  is  important  to  note  that  the  detector  occludes  the  source 
when  it  is  near  the  backscattering  position.  This  results  in  a  loss  of  3°  to  6°  of  scattering 
amplitudes  centered  about  the  source  position.  These  data  points  are  not  plotted. 

The  final  desired  measurement  is  the  mean  scattered  pressure  as  a  function  of  scattering 
angle,  given  a  fixed  angle  of  incidence.  Since  the  scattering  mechanisms  working  at  the 
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interface  are  strongly  a  function  of  frequency,  we  are  interested  in  measuring  the  scattered 
field  at  specific  frequencies.  Two  different  approaches  to  this  problem  can  be  taken.  One  is  to 
record  the  scattered  energy  given  a  broad  band  incident  wavelet  and  then  to  decompose  the 
recorded  energy  as  function  of  frequency.  However,  since  it  is  unclear  exactly  what  numerical 
effects  may  be  introduced  with  the  narrow  oand  filtering  we  take  the  second  approach  which 
involves  no  filtering  and  gives  cleaner  results.  We  introduce  directly  via  the  transducer  source 
a  monochromatic  wave  of  given  frequency.  The  continuous  sinusoidal  wave  is  approximated 
well  by  a  finite  sinusoid  ranging  from  35  to  100  cycles.  The  final  constant  amplitude  of  the 
scattered  pressure  is  then  recorded. 

A  single  realization  of  scattering  from  the  irregular  interface  results  in  large  fluctuations 
in  scattered  pressure  as  a  function  of  the  scattering  angle.  It  is  necessary,  therefore,  to 
average  experimentally  as  we  averaged  numerically,  so  as  to  obtain  a  final  mean  reflection 
coefficient.  In  optical  theory,  averaging  is  accomplished  by  illuminating  a  field  lens  that 
is  much  larger  than  these  fluctuations.  The  fluctuations,  referred  to  as  speckle  in  optical 
terminology  (O’Donnell  and  Mendez,  1987),  are  integrated  over  a  specified  solid  angle  giving 
the  average  intensity  scattered  at  that  angle.  In  seismic  experiments  it  is  very  difficult  to  use 
an  integrating  lens.  This  difficulty  arises  mainly  from  the  limited  size  and  frequency  range 
of  ultrasonic  experimentation.  As  a  result,  we  chose  to  follow  the  numerical  approach  and 
to  create  different  independent  realizations  of  the  interface.  The  pressure  of  waves  scattered 
from  an  ensemble  of  interface  realizations  was  averaged  to  determine  the  final  mean  scattered 
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pressure  at  each  scattering  angle.  The  computer-controlled  arm  allowed  for  rapid  acquisition 
ot  data  with  the  reproducibility  required  for  this  averaging  scheme.  A  typical  measurement 
consisted  of  measuring  the  scattered  pressure  from  an  independent  statistical  realization  of 
the  interface.  Each  independent  realization  was  acquired  by  rotating  and  shifting  the  sample 
in  the  sequence  shown  in  Figure  8.  Since  rotating  the  surface  with  respect  to  the  incident 
be  n  iormed  another  scattering  geometry,  many  different  realizations  of  the  interface  were 
obtained.  In  general,  the  sample  was  rotated  by  60°  staggered  steps,  followed  by  1.25  cm 
shifts  of  the  block,  which  placed  the  axis  of  rotation  for  the  source  and  receiver  on  concentric 
circles  about  the  center  of  the  glass  cylinder.  This  movement  sequence  can  easily  give  more 
than  150  different  realizations  of  the  interface. 

Due  to  interest  in  the  enhancement  of  energy  travelling  directly  back  .'ovvards  the  source, 
it  was  desirable  to  create  a  source-receiver  design  that  retrieved  energy  in  the  occluded  zone 
near  the  source.  This  was  achieved  by  constructing  a  four-layered  piezo-film  detector.  The 
general  idea  behind  the  piezo-fdm  receiver  is  straightforward.  A  single  layer  of  film  has  an 
impedance  very  similar  to  water,  and,  as  a  result,  when  the  piezo-film  is  placed  directly  in 
front  of  the  source,  the  source  energy  transmits  almost  completely  through  the  detector, 
allowing  most  of  the  energy  to  travel  towards  the  interface.  The  incident  wave  is  then 
diffracted  from  the  interface,  and  the  scattered  energy  travelling  directly  back  towards  the 
source  is  recorded  as  it  transmits  a  second  time  through  the  piezo-film  receiver.  A  four- 
layer  piezo-film  receiver  was  constructed  using  conducting  glue  and  in-parallel  connections, 


49 


significantly  increasing  the  signal  to  noise  ratio.  The  construction  of  this  receiver  is  discussed 
in  detail  in  Appendix  B.  It  is  important  to  note  that  the  piezo-film  receiver  is  sensitive  to 
energy  arriving  both  from  in  front  and  from  behind  the  receiver.  As  a  result,  extreme  care 
was  taken  to  choose  a  time  window  of  analysis  that  did  not  contain  multiple  scattering  from 
the  source  and  the  tank  wall.  In  addition,  as  we  will  show,  piezo-film  of  this  thickness  does 
have  a  substantial  reflection  coefficient,  therefore,  a  resonance  between  source  and  receiver 
had  to  be  avoided,  too. 

We  stress  that  the  data  discussed  in  the  next  section  represents  the  mean  diffusely  scat¬ 
tered  signal,  as  a  function  of  angle,  for  a  fixed  solid  angle  of  data  acquisition.  No  artificial 
angular  factors  are  introduced  to  the  data  even  though  the  apparent  vertical  wavelength, 
acting  at  the  surface  in  this  experiment,  varies  inversely  with  the  cosine  of  the  scattering 
angle.  We  do  not  place  any  absolute  vertical  scales  on  the  data,  although  for  comparison 
the  data  are  normalized  to  the  numerical  RC  calculated  at  normal  incidence. 


Surface  Scattering  Measurements 

In  this  section  we  present  the  average  reflection  coefficient  measurements  obtained  using  the 
roughest  glass  surface  fabricated.  This  surface  was  chosen  because  many  scattering  mecha¬ 
nisms  will  play  their  strongest  role.  Figure  9a  shows  the  target  surface  height  distribution, 
independent  of  lateral  position,  and  Figure  9b  shows  the  histogram  of  the  surface  height, 
based  on  surface  profilometer  measurements  of  the  actual  surface.  The  histogram  shows 
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that  the  data  nicely  matches  a  Gaussian  probability  distribution  with  a  standard  deviation 
of  approximately  0.6  mm.  This  is  close  to  the  target  value,  8  =  0.71.  Figure  10a  gives  the 
target  Gaussian  autocorrelation  function,  and  Figure  10b  shows  the  actual  autocorrelation 
function  calculated  from  profilometer  measurements.  As  previously  predicted,  the  correla¬ 
tion  length  of  1.4  mm  is  greater  than  the  target  value  of  1  mm.  Also  plotted  are  Gaussian 
and  exponential  autocorrelation  functions  with  the  same  correlation  length  as  the  data.  The 
autocorrelation  function  is  very  close  to  a  Gaussian  correlation  function  at  the  more  impor¬ 
tant.  smaller  lags.  At  larger  lag  distances,  the  surface  lies  directly  in  between  a  Gaussian 
and  an  exponential  correlation  function. 

Figure  11a  gives  a  3-D  grayscale  plot  of  the  irregular  surface  based  on  the  surface  pro¬ 
filometer  measurements.  Figure  lib  plots  the  surface  height  for  a  profile  taken  across  the 
surface,  while  Figure  11c  shows  a  numerically  generated  Gaussian  and  exponential  surface 
given  the  same  standard  deviation  and  correlation  length.  It  is  clear  that  the  Gaussian 
surface  matches  the  experimental  interface  well,  both  in  the  observed  slopes  and  the  lateral 
scale  of  the  irregularities.  Since  the  frequency  domain  representation  of  the  surface  goes 
as  the  square  root  of  its  correlation  function,  the  exponential  surface  should  contain  larger 
amounts  of  energy  at  both  low  and  high  frequencies.  This  is  seen  clearly  since  the  exponen¬ 
tial  surface  contains  lower  amplitude  short  wavelength  irregularities  that  were  not  observed 
on  the  experimental  surface.  Although  the  Gaussian  gives  a  good  fit  to  the  experimental 
data,  out  of  interest,  we  shall  still  plot  the  results  for  an  exponential  surface. 
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Based  on  the  measurements  above,  the  final  glass  interface  slopes  steeply  at  approxi¬ 
mately  30°  rms,  and  the  impedance  contrast  at  the  fluid-glass  interface  is  large  as  the  glass 
interface  has  properties  very  similar  to  an  igneous  material  (see  Figure  3).  As  a  result, 
multiple-scattering  and  shadowing  effects  can  play  a  significant  role  at  both  small  and  large 
incident  angles,  and  approximate  linear  theories,  such  as  the  KirchhofF  and  Born  approaches, 
break  down.  Therefore,  the  mean  scattered  pressure  measured  experimentally  is  compared 
with  the  reflection  coefficients  calculated  with  the  Boundary  Integral  technique,  formulated 
earlier.  This  technique  includes  all  multiple  scattering  and  wave  conversions  at  the  interface. 

Case:  A  =  0.71a 

Figure  12  shows  one  realization  of  the  interface  given  an  incident  pulse  with  a  center  frequency 
of  1.5  MHz  and  a  half-power  width  of  0.25MHz.  This  realization  corresponds  to  a  beam 
impinging  on  the  surface  with  an  incident  angle  of  20°.  The  source  pulse  is  shown.  This 
pulse  is  alo  shown  reflected  from  a  plane  interface,  in  which  case  energy  arrives  only  in 
the  specular  direction.  Referring  to  the  polar  seismogram,  it  is  clear  that  the  large  scale 
surface  roughness  scatters  energy  over  most  forward  and  back  scattering  angles.  The  energy 
is  spread  over  a  large  time  interval  and  amplitudes  vary  rapidly  as  a  function  of  scattering 
angle.  In  general,  it  is  difficult,  given  this  single  model,  to  determine  quantitatively  which 
scattering  mechanisms  are  working  at  the  surface. 

Our  first  continuous  wave  analysis  is  carried  out  at  1.5  MHz,  the  center  frequency  of 
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the  seismogram  above,  so  that  A  =  .71a  =  1.00  mm.  Figure  19  shows  one  experimental 
realization  of  the  fluid-glass  surface  at  each  of  four  incident  beam  angles:  0°,  20°,  30°,  and 
60°.  Clearly  the  amplitudes  in  each  realization  vary  strongly  as  a  function  of  scattering 
angle,  6S.  As  discussed  earlier,  these  fluctuations  can  be  removed  by  averaging  over  a  finite 
number  of  realizations  and  obtaining  a  mean  reflection  coefficient.  The  total  mean  reflection 
coefficients  for  both  a  Gaussian  and  an  exponential  surface  are  given  in  Figures  14-17  and 
Figures  20-23.  At  the  bottom  of  these  figures,  the  experimental  mean  reflection  coefficients 
are  given  along  with  the  SD  of  the  finite  average,  showing  the  deviation  of  these  reflection 
coefficients  from  the  final  mean  reflection  coefficient  corresponding  to  a  full  ensemble  of 
realizations.  Negative  scattering  angles  (6S  <  0)  correspond  to  backscattering  in  all  plots. 
We  also  stress  that  given  the  incident  wavelength,  the  surface  is  extremely  irregular  and  the 
specular  reflection  is  largely  disrupted. 

Figure  14  shows  the  total  mean  scattered  pressure  as  a  function  of  scattering  angle  given 
a  normally  incident  acoustic  beam.  Upon  comparing  the  numerical  and  experimental  data  it 
is  clear  that  the  2-D  numerical  results  for  a  Gaussian  interface  match  the  3-D  scattered  data 
extremely  well.  The  fluctuations  in  the  data  are  mostly  within  one  standard  deviation  of  the 
finite  average.  There  are  a  number  of  interesting  aspects  in  the  curves.  The  experimental 
data  shows  a  strong  peak  amplitude  propagating  back  towards  the  source  at  0S  =  0°;  this 
is  predicted  well  by  the  numerical  reflection  coefficient.  The  width  of  this  peak  is  approxi¬ 
mately  35°.  Notice  that  there  is  considerable  scattering  at  all  angles,  dropping  off  linearly 


with  increasing  scattering  angle.  The  exponential  surface  also  predicts  the  general  form  of 
scattering  well,  but  the  higher  frequency  component  appears  to  destroy  the  enhancement  of 
energy  scattered  back  towards  the  source. 

A  similar  form  of  scattering  is  exhibited  in  Figures  15  and  16,  which  show  the  mean 
scattered  pressure  for  an  incident  angle  of  20°  and  30°,  respectively.  Both  the  numerical 
results  for  the  Gaussian  interface  and  the  experimental  data  match  very  well.  There  are  also 
two  remarkable  aspects  to  these  curves.  First,  on  both  reflection  coefficients  two  peaks  can  be 
identified  by  eye.  One  is  a  broad  peak  occurring  in  the  forward  scattered  direction.  The  other 
is  much  narrower  and  occurs  in  the  retroreflective  direction  0,  =  —6q.  This  “retroreflective” 
peak  loses  amplitude  as  the  incident  angle  is  increased,  sinking  further  into  the  surrounding 
reflection  coefficient  curve.  This  reUoreflectance  is  clearly  supported  by  the  ultrasonic  data. 
Second,  both  curves  become  strongly  asymmetric.  However,  upon  comparing  the  curves,  the 
2-D  numerical  model  shows  more  backscattering  and  less  forward  scattering  than  the  3-D 
ultrasonic  data.  This  trend  becomes  more  prominent  as  the  incident  angle  is  increased.  The 
exponential  curve  also  follows  the  experimental  reflection  coefficient  closely,  showing  again 
a  less  distinct  retroreflective  peak. 

Figure  17  gives  the  mean  scattered  pressure  for  a  beam  incident  at  60°.  In  this  case  there 
are  no  distinct  signs  of  enhanced  backscattering  in  either  the  numerical  or  the  experimental 
data.  However,  energy  is  scattered  uniformly  over  most  backscattering  angles.  This  energy 
does  not  drop  off  until  the  retroreflective  angle  is  exceeded  in  the  backscattering  region. 
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Most  dramatic  is  the  continuation  of  the  trend  observed  at  the  smaller  incident  angles  above. 
Specifically,  the  numerical  data  clearly  shows  more  backscattering  and  less  forward  scattering 
than  the  experimental  data. 

Case:  A  =  2.14a 

Figure  18  shows  one  realization  of  the  interface  given  an  incident  pulse  with  a  center  frequency 
of  0.5  MHz  and  a  half-power  width  of  250  kHz.  This  realization  corresponds  to  a  beam 
impinging  on  the  surface  at  a  20°  incident  angle.  The  pulse  reflected  from  a  plane  interface 
is  shown,  with  the  energy  again  arriving  only  in  the  specular  direction.  Referring  to  the 
scattered  seismogram,  it  is  clear  that,  even  at  this  lower  frequency,  energy  is  scattered  over 
most  forward  and  backscattered  angles.  Given  this  one  deterministic  case,  it  is  difficult  to 
determine  quantitatively  the  scattering  mechanisms  operating  at  this  frequency,  or  to  define 
how  they  might  differ  from  those  in  the  higher  frequency  case. 

This  second  continuous  wave  analysis  was  carried  out  at  0.5  MHz,  the  center  frequency 
of  the  seismogram  above,  so  that  A  =  2.14a  =  3.0  mm.  Figure  13  shows  one  experimental 
realization  of  the  fluid-glass  surface  at  each  of  the  four  incident  beam  angles.  Once  again 
the  amplitudes  for  each  realization  vary  strongly  as  a  function  of  scattering  angle,  although 
not  as  strongly  as  the  A  =  .71a  case.  Figure  20  shows  the  comparison  between  the  averaged 
numerical  and  experimental  reflection  coefficients  given  a  normally  incident  beam.  The  2-D 
numerical  results  for  a  Gaussian  interface  predict  the  experimental  observations  well.  All  of 
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the  experimental  data  sits  within  the  standard  deviation  of  the  finite-average.  Comparing 
these  curves  to  the  curves  for  A  =  0.71a,  a  number  of  distinct  differences  are  apparent. 
Most  noticeable  is  the  widening  of  the  retroreflective  peak  width  from  about  35°  to  greater 
than  60°.  This  widening  is  apparent  in  both  the  experimental  and  the  numerical  data.  The 
reflection  coefficient  for  an  exponential  interface  again  shows  much  lower  retroreflectance 
than  for  the  Gaussian  interface. 

Figures  21  and  22  both  show  that  the  numerical  results  over  a  Gaussian  interface  predict 
the  asymmetric  trends  in  the  experimental  data  very  well  for  incident  angles  of  20°  and 
30°,  respectively.  However,  distinct  differences  do  occur  between  the  two  curves.  First,  as 
the  incident  angle  increases,  the  2-D  numerical  results  again  show  more  backscattering  and 
less  forward  scattering  than  the  3-D  ultrasonic  data.  A  broad  retroreflective  peak  is  both 
predicted  and  observed  at  20°  and  30°  incidence,  supporting  the  existence  of  retroreflectance. 
Unfortunately,  the  height  of  these  peaks  are  of  the  same  order  as  the  standard  deviation  of  the 
experimental  average,  not  allowing  for  a  direct  verification  of  retroreflectance.  Numerically, 
the  exponential  interface  does  give  rise  to  a  retroreflective  peak;  this  peak  is  smaller  than 
the  peak  predicted  by  the  Gaussian  surface. 

Figure  23  shows  the  mean  reflection  coefficient  for  an  incident  angle  of  60°.  In  this  case 
enhanced  backscattering  is  not  predicted  numerically  or  observed  experimentally.  At  this 
lower  frequency,  the  2-D  numerical  model  predicts  more  backscattering  and  less  forward 
scattering  than  the  3-D  ultrasonic  data.  In  addition,  the  numerical  model  predicts  a  much 
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smaller  specular  peak  than  is  observed  experimentally.  Although  the  amplitudes  are  different, 
the  numerical  curves  do  predict  well  the  uniform  scattering  of  energy  into  the  fluid  above  as 
shown  in  the  data. 

Retroreflectance  is  clearly  supported  by  the  ultrasonic  data  above.  However,  the  retrore- 
flective  peak  height  is  still  on  the  same  order  as  the  standard  deviation  of  the  finite-average 
in  each  case.  This  makes  it  difficult  to  verify  the  existence  of  “enhanced  backscattering.”  For 
this  reason,  data  was  recorded  near  the  retroreflective  direction  using  the  partially  transpar¬ 
ent  piezo-film  receiver.  The  experimental  procedure  is  described  in  Appendix  B,  the  steps 
obtaining  the  final  reflection  coefficient  are  summarized  in  Figure  24,  and  the  final  average 
RC  observed  with  the  piezo-film  receiver  is  superposed  on  Figure  21.  The  data  has  been 
scaled  to  the  amplitudes  received  with  the  flat-bottomed  transducers.  The  scattered  pres¬ 
sure  was  measured  between  the  backscattering  angles  of  40°  and  5°  (—40°  <  9S  <  —5°), 
and  65  surface  realizations  were  averaged.  In  this  case,  the  size  of  the  retroreflective  peak  is 
larger  than  the  corresponding  SD  of  the  average.  The  2-D  numerical  model  predicts  the  3-D 
experimental  data  very  well,  to  within  the  standard  deviation  of  the  finite  average.  A  dis¬ 
tinct  peak  is  observed  in  the  retroreflective  direction  with  a  slightly  narrower  form  than  the 
numerically  generated  peak.  This  result  strongly  supports  the  enhancement  of  backscattered 
energy  due  to  multiple  scattering  from  the  glass  interface. 
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General  Discussion 


Generally,  the  numerically  derived  mean  reflection  coefficients  calculated  over  an  acoustic- 
elastic  interface  show  retroreflective  trends  similar  to  those  observed  for  the  SH  and  P-SV 
cases.  First,  the  width  that  the  retroreflective  peak  appears  to  be  is  consistent  with  the 
multiple  scattered  constructive  phase  argument  summarized  in  the  introduction.  In  this 
case,  when  the  wavelength  is  increased  by  a  factor  of  three,  both  the  numerically  derived 
curve  and  the  experimental  data  show  a  factor  of  three  increase  in  peak  width,  from  35°  at 
A  =  ,71a  to  greater  than  60°  at  A  =  2.14a.  Second,  as  the  incident  angle  is  increased,  the 
retroreflective  peak  amplitude  tends  to  decrease  relative  to  the  remaining  portion  of  the  DRC. 
The  retroreflective  peak  in  both  the  experimentally  and  numerically  derived  curves  seems  to 
disappear  around  an  incident  angle  equal  to  the  30°  rms  slope  of  the  interface.  Although 
not  studied  directly  here  it  seems  likely,  based  on  the  work  of  Schultz  (1993a),  that  the 
retroreflective  peak  height  will  tend  to  decrease  as  the  impedance  contrast  is  lowered  and 
more  energy  is  allowed  to  penetrate  the  interface.  Along  the  same  lines,  the  retroreflective 
peak  amplitude  is  likely  to  diminish  as  the  rms  slope  of  the  interface  is  decreased,  since  not 
as  many  time-reversed  paths  can  be  obtained  with  lower  slopes. 

Interestingly,  as  the  incident  angle  was  increased,  the  numerical  results  above  consistently 
predicted  “more  backscattering  and  less  forward  scattering”  than  observed  in  the  3-D  exper¬ 
imental  results.  We  stress  that  the  amplitudes  for  the  2-D  and  3-D  are  normalized  to  each 
other  at  normal  incidence,  so  that  the  absolute  amplitude  of  backscattering  is  not  given.  The 
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above  statement,  “more  backscattering  and  less  forward  scattering,”  means  that  backscatter- 
ing  decreases  and  forward  scattered  amplitudes  increase  more  slowly  in  the  2-D  case  than  in 
the  3-D  case  as  the  incident  angle  is  increased.  Two  explanations  may  clarify  this  deviation. 
First,  the  glass  interface  used  in  the  experiment  may  not  follow  Gaussian  statistics  exactly. 
Therefore,  the  observed  differences  may  be  a  direct  result  of  differences  in  the  statistics  of 
the  model.  However,  the  interface  statistics  were  well-constrained  using  surface  profilometer 
measurements.  A  second  possible  explanation,  is  that  there  may  be  distinct  differences  in 
the  scattering  due  to  an  inherent  difference  between  2-D  and  3-D  scattering  mechanisms.  In 
the  case  of  2-D  scattering,  a  peak  along  the  interface  is  actually  an  infinite  ridge  extending 
in  the  x2-direction.  A  wave  incident  on  the  side  of  this  ridge  has  only  three  probable  routes 
of  getting  to  the  receiver  located  on  the  opposite  side  of  the  ridge.  The  wave  can  either 
multiply  scatter  within  the  valley,  transmit  through  the  ridge,  or  diffract  over  the  very  peak 
of  the  ridge.  Each  of  these  paths  exists,  but  a  wave  loses  a  large  amount  of  energy  along  any 
of  these  paths.  In  the  3-D  case,  the  surface  has  one  more  degree  of  freedom  so  that  a  peak 
along  the  interface  can  vary  in  all  directions.  Therefore,  at  normal  incidence,  out  of  plane 
scattering  allows  energy  to  arrive  randomly  from  all  directions  and  may  increase  the  amount 
of  observed  backscattering.  As  the  incident  angle  is  increased,  energy  travels  also  out  of  the 
incident  plane.  However,  unlike  the  2-D  interface,  the  3-D  nature  of  the  interface  may  allow 
energy  to  pass  around  obstructing  peaks,  reducing  the  amount  of  backscattered  energy  in 
the  measure  that  the  incident  angle  is  increased.  For  instance,  energy  may  diffract  from  the 
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flank  of  a  nearby  peak  and  propagate  in  the  forward  direction  back  into  the  receiver  plane. 
In  this  case,  out-of-plane  scattering  would  work  to  reduce  the  amount  of  backscattering,  and 
additional  energy  may  be  allowed  to  propagate  as  forward  scattered  energy,  explaining  the 
difference  between  the  2-D  and  3-D  trends. 

Conclusions 

In  this  study,  we  were  able  to  generate  within  reasonable  accuracy  a  3-D  characterized  ran¬ 
dom  interface  with  Gaussian  statistics.  An  interface  with  approximately  a  Gaussian  surface 
height  distribution  and  Gaussian  correlation  function  was  generated  using  a  glass  etching 
procedure  and  photoresist  templates.  The  resulting  surface  distribution  was  confirmed  using 
surface  profilometer  measurements.  Scattered  pressures  were  then  acquired  over  this  surface 
and  compared  directly  to  numerical  results  calculated  over  a  2-D  interface  with  the  same 
statistical  parameters. 

Specifically,  we  have  shown  that  2-D  numerical  models  of  an  acoustic-elastic  interface 
with  Gaussian  statistics  predict  the  3-D  scattering  that  is  observed  experimentally  very 
well,  when  the  incident  wave  is  near  normal  incidence.  Numerical  results  predict  the  large 
amount  of  observed  incoherent  backscattering  and  forward  scattering.  Numerical  results 
also  predict  the  elimination  of  the  specular  reflection  at  many  incident  angles.  Most  striking 
is  the  prediction  of  a  "retroreflective”  amplitude  peak.  The  numerical  results  predict  that 
the  peak  amplitude  decreases  as  the  incident  angle  increases  and  the  peak’s  width  is  directly 
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proportional  to  the  ratio  of  the  incident  wavelength  and  the  correlation  length  of  the  interface. 
Experimentally,  enhanced  backscattering  is  strongly  supported  by  observations  at  normal 
incidence  and  at  20°-incidence.  The  peak’s  amplitude  appears  to  decrease  dramatically  when 
the  incident  angle  becomes  greater  than  the  rms  slope  of  the  interface.  Numerical  modeling 
of  an  exponential  surface  with  the  same  correlation  length  and  standard  deviation  as  the 
glass  interface  gives  results  very  similar  to  the  results  of  the  Gaussian  surface  although  the 
results  show  more  general  backscattering.  Retroreflectance  from  the  exponential  interface 
was  difficult  to  identify  in  almost  all  cases  studied. 

The  2-D  numerical  curves  deviate  from  the  experimental  curves  at  incident  angles  greater 
than  30°.  Specifically,  as  the  incident  angle  is  increased,  the  2-D  models  predict  more 
backscattering  and  less  forward  scattering  than  observed  in  the  experimental  data.  This 
trend  becomes  much  stronger  as  the  incident  angle  approaches  grazing  angles.  As  the  surface 
in  this  experiment  is  well  characterized,  this  appears  to  result  from  a  distinct  difference 
between  2-D  and  3-D  scattering  mechanisms.  In  this  case  the  peaks  and  valleys  of  3-D 
surfaces  appear  to  backscatter  less  energy  into  the  plane  of  incidence  than  the  ridges  along 
2-D  interfaces  when  the  incident  angle  is  increased.  In  addition,  they  allow  more  energy  to 
scatter  forward  into  the  plane  of  incidence.  Even  though  the  amplitudes  did  not  match  at 
larger  incident  angles,  the  trends  of  the  reflections  coefficients  did  match.  Both  numerical 
and  experimental  results  showed  energy  scattered  uniformly  over  most  scattering  angles  as 
a  60°  incident  angle  was  approached.  This  results  in  negative  phase  velocities,  large  phase 
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velocities,  and  a  large  amount  of  interference  in  seismic  data  recorded  in  the  fluid  above  the 
interface. 

Although  tests  were  not  carried  out  on  epoxy  surfaces  in  this  study,  an  epoxy  surface 
was  generated  for  the  purpose  of  profilometer  measurements.  Essentially,  the  irregular  glass 
surface  was  used  as  a  mold.  After  adding  a  separating  solution,  an  epoxy  model  was  generated 
and  separated  from  the  glass  mold.  We  have  concluded  that  the  glass  surface  can  be  used  to 
make  models  out  of  different  substances,  each  with  different  material  properties.  Future  work 
will  therefore  include  studies  to  determine  how  the  material  properties  affect  scattering  from 
surfaces  with  identical  height  distributions.  Although  in  this  paper  we  chose  to  fabricate  a 
Gaussian  surface,  the  etching  process  used  here  can  also  be  used  to  create  an  interace  with 
exponential  statistics. 
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Appendix  A:  Numerical  Formulation — Acoustic-Elastic  Case 

We  numerically  solve  the  integral  equations  (20)  and  (21)  by  first  integrating  them  over  a 
finite  interval,  L,  then  converting  them  to  a  set  of  3 N  linear  equations.  Separating  each  of 
the  integrals  into  a  sum  of  N  integrals,  each  integrated  over  an  increment  Ax,  centered  at 
the  interface  points  we  can  express  (21),  which  represents  the  scattered  displacement  in  the 
lower  medium,  as  a  sum  over  integrals  centered  at  the  points 

xn  =  -L/2  +  (n  —  ^)Ax,  n  =  1,2,3, ...,  N.  ( A.l) 

Evaluating  the  total  displacement  at  the  center  of  each  element,  xm,  the  integral  approxi¬ 
mations  can  be  written  as 

(xmlx[),  (A.2) 

2 

r  Xrx  “f*  --1- 

D$){xm\xn)  =  /  J  dx\D?\xm\x\), 

''Xrx - 2 

which  we  solve  for  kk\x  very  small. 

In  the  acoustic  medium  the  integrals  which  we  must  approximate  analytically  can,  as¬ 
suming  that  the  source  amplitude  functions  24  are  slowly  varying  along  the  interface,  be 
written  as 

rxn+ 

r‘/)(.rm|xn)  =  /  J  dx\Pf)(xm |xi),  (A. 3) 

''Xrx  2 

D\i\xm |xn)  =  f  dx1Di/\xm\x'1), 

X  "2  "' 
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and  it  can  be  easily  shown  by  letting  x3  =  C(xi)  +  e  and  integrating  all  terms  corresponding 
to  first  order  in  k^ Ax  in  the  limit  as  t  — ♦  0+  that  the  boundary  integral  formulation  in 
the  fluid  can  be  represented  as  a  simple  combination  of  both  the  elastic  approximation  and 
the  SH  approximation  (Schultz  and  Toksoz,  1993a).  In  the  acoustic  medium  the  Green’s 
function  terms  are  written  as 

D(/\x\x')  =  A  </>4/)2G(/)(x;x>,-  (A.4) 

r(/)(x  |x')  =  T{SH){x \x')[h->\] 

where  the  traction  term  is  identical  to  the  SH  traction  term  with  A^  in  the  fluid  is  substituted 
for  /z  in  the  SH  case.  In  the  elastic  medium  these  terms  can  be  written  as 

0<,"(!sls')  =  /‘(,|CS(s;s>,  (A.5) 

r"(,)(  sis')  =  r"(‘’(s|s') 

where  the  traction  term  is  just  that  of  the  P-SV  case.  The  displacement  term  is  a  function  of 
the  elastic  Green's  function  developed  in  the  P-SV  formulations.  In  the  case  of  integrating 
over  singularities,  expressing  the  integrands  as  a  combination  of  Taylor  and  asymptotic  series 
expansions,  keeping  terms  up  to  order  kj^  Ax,  and  integrating  over  Ax  gives  a  solution  which 
reduces  to  a  simple  combination  of  the  SH  and  P-SV  solutions.  This  includes  the  first  term 
in  A.5  which  multiplies  an  additional  normal  function.  A  direct  substitution  from  Schultz 
and  Toksoz  ( 1993a)  and  Schultz  and  'oksoz  (1993b)  gives  the  final  linear  system  of  equations 
which  is  then  solved. 
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Appendix  B:  Ultrasonic  Design 


All  experiments  were  done  in  our  in-house  ultrasonic  water  tank  laboratory.  Ultrasonic 
water  tank  modelling  is  a  powerful  tool  as  real  field  situations  can  be  scaled  down  by  four  to 
five  orders  of  magnitude  and  studied  in  a  well  controlled  laboratory  environment.  The  first 
section  of  this  appendix  describes  the  ultrasonic  configuration.  The  second  section  describes 
in  detail  the  construction  of  a  four  layer  piezo-film  receiver. 

Ultrasonic  Water  Tank  Design 

A  diagram  of  the  experimental  setup  is  shown  in  Figure  B.l.  A  solid  elastic  model  is 
submerged  in  a  water  tank  measuring  1.0  m  by  0.5  m  by  0.5  m  in  height.  An  ultrasonic  wave 
is  generated  at  a  given  source  transducer  with  an  input  voltage  from  a  Hewlett  Packard 
1048 A  function  generator.  For  voltage  output  above  +/  —  5  V  a  Hewlett  Packard  467A 
Power  Amplifier  was  utilized  to  obtain  a  +/  —  10  V  of  output.  The  source  function  generator 
was  used  to  create  a  finite  length  sine  wave  with  a  specific  frequency  and  a  finite  number 
of  cycles.  The  input  wave,  received  with  a  corresponding  piezoelectric  material,  is  amplified 
from  40  to  60  db  with  a  Panametrics  5660B  preamplifier. 

In  the  case  of  the  finite  pulse  seismograms,  a  Krohn-Hite  3202R  high-low  cutoff  filter  was 
applied  to  the  input  signal.  In  the  continuous  wave  study  the  input  voltage  was  increased  and 
the  filter  bypassed  to  remove  any  distortion  it  may  cause.  Signal  to  noise  ratios  were  improved 
by  stacking  recorded  signals  (usually  8  stacked  shots).  The  final  signal  was  then  digitized  by 
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a  Data  Precision  DATA  6000  digital  oscilloscope  with  12-bit  amplitude  resolution.  Although 
the  system  can  simulate  a  full  3-D  seismic  study  using  six  automated  step  motors,  this  study 
was  carried  out  with  only  one  motor.  The  receiver  was  connected  rigidly  to  a  rotating  arm 
and  then  rotated  in  a  semicircle  about  a  rotation  axis  in  .9°  steps.  The  digitizer  and  the  step 
motor  controller  are  interfaced  with  an  IBM  PC- AT  computer  through  an  IEEE-488  interface 
bus.  The  final  digitized  data  is  stored  in  the  IBM  PC-AT  and  following  the  experiment,  is 
transferred  to  a  Digital  DECstation  5000/25  machine.  The  total  time  required  to  record  and 
transfer  the  pressure  recorded  for  a  single  realization  of  an  interface  was  on  the  order  of  one 
hour. 

Piezo-Film  Receiver 

In  designing  an  experiment  which  can  study  the  energy  scattered  at  and  near  retroreflective 
angles  (energy  travelling  back  towards  the  source),  one  must  be  able  to  place  a  receiver  very 
close  to  the  source.  There  are  a  number  of  approaches  which  allow  recording  directly  in  this 
region.  One  approach  is  to  make  the  source  as  small  as  possible  so  that  when  the  receiver 
gets  near  the  source  it  does  not  interfere  with  the  incident  beam.  The  smaller  the  source 
and  receiver  the  more  the  information  recorded  at  retroreflective  angles.  A  second  approach 
is  to  use  a  focusing  beam  which  has  its  focus  at  the  receiver.  In  this  case  the  beam  is  its 
narrowist  at  the  receiver  allowing  recording  very  near  the  source  angle.  A  third  approach  is 
to  use  an  acoustic  beam  splitter,  which  reflects  incident  energy,  coming  from  the  side  of  the 
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experiment,  towards  the  interface  and  then  allows  energy  to  pass  back  through  the  splitter 
on  its  second  pass  to  the  receiver.  The  fourth  approach  is  to  create  a  receiver  which  is 
transparent  to  acoustic  energy  giving  one  free  reign  to  record  about  the  source  transducer. 
Unfortunately,  the  first  two  approaches  are  not  applicable  to  this  study.  The  first  option  is 
not  practical  since  the  size  of  the  source  transducer  is  required  to  be  larger  than  a  certain 
limit  to  propagate  a  narrow  acoustic  beam.  The  second  option  is  also  difficult  to  implement 
because  of  the  water  tank  dimensions.  If  the  beam  is  focused  on  a  receiver  located  only  a 
short  distance  away  the  beam  spreading  becomes  too  large. 

YVe  therefore  chose  the  fourth  approach.  We  describe  a  piezo-film  transducer  that  is 
essentially  transparent  to  acoustic  energy.  Figure  B.2(a,c)  shows  a  detailed  view  of  the 
piezo-film  receiver  configuration.  In  the  case  of  a  single  piezo-film  layer,  the  film  is  placed 
between  two  silver  conducting  layers.  Since  the  signal  received  comes  only  from  areas  covered 
by  conducting  metal,  it  is  advantageous  to  apply  the  metal  coating  as  a  liquid  glue  which 
allows  for  full  control  over  the  receivers  shape  and  its  sensitivity.  In  this  study  the  conducting 
element  was  given  the  same  circular  shape  as  the  flat-bottomed  transducers  described  earlier. 
Due  to  the  low  sensitivity  of  one  layer  of  piezo-film,  it  is  desirable  to  stack  multiple  layers. 
However,  t  here  is  a  balance  between  the  thickness  of  the  receiver  and  its  impedance  contrast. 
As  more  piezo-film  layers  are  stacked  the  impedance  of  the  stack  increases  quickly.  For  this 
study,  a  four  layer  stack  shown  in  Figure  B.2(c)  using  110/i  film  was  optimal.  The  final 
thickness  of  the  piezo  film  is  approximately  550/x  thick  with  a  transmis.:  .efficient  at 
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normal  incidence  of  approximately  0.7S.  This  leaves  the  piezo-film  with  a  .22  reflection 
coefficient.  As  a  result,  reflected  energy  does  reverberate  between  the  source  and  receiver 
retaining  a  very  high  amplitude.  A  recording  window  was  chosen  such  that  it  did  not  include 
any  of  this  reverberative  energy. 

The  final  sheet,  measuring  .25  m  by  .12  m,  was  supported  around  the  edges  using  a  thin 
frame  (Figure  B.3).  In  addition  to  supporting  the  film,  the  frame  was  used  to  bend  the 
sheet  so  it  followed  the  curvature  of  the  semicircle  about  which  the  receiver  is  stepped.  This 
guaranteed  that  the  source  beam  was  normally  incident  on  the  film  during  receiver  rotation, 
therefore  reducing  any  source  beam  distortion  created  by  oblique  incidence. 

As  this  is  a  new  transducer  design,  the  receiver  properties  were  studied.  Figure  B.4 
shows  the  receiver  sensitivity  pattern  as  a  function  of  incident  wave  angle,  with  an  incident 
frequency  of  500  kHz.  The  stacked  piezo- film  receb  '  as  a  very  directed  sensitivity  pattern 
which  is  similar  in  nature  to  the  flat-bottomed  i  acers  with  the  same  diameter.  The 
piezo-film  also  appears  to  have  a  broad  frequency  response  making  it  a  very  powerful  tool.  In 
addition,  this  element  may  be  used  successfully  in  rock  physics  to  record  on  curved  surfaces 
by  gluing  the  film  directly  to  the  rock  surface.  This  avoids  the  problem  of  transducer  coupling 
encountered  with  a  flat-bottomed  transducer. 
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Figure  1:  Peak- valley  sequence  along  an  interface  showing  an  example  of  the  time-reversed 
paths  which  lead  to  enhanced  backscattering.  The  solid  line  shows  a  forward  scattered 
path,  while  the  dashed  line  shows  the  corresponding  time-reversed  path.  These  two  paths 
constructively  interfere  to  give  an  increase  in  amplitude  diffracted  back  in  the  direction 
of  the  source. 
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Figure  2:  The  geometry  used  to  formulate  the  numerical  model  of  scattering  from  a  highly 
irregular  interface.  The  upper  acoustic  volume  is  separated  from  the  lower  elastic  medium 
by  a  highly  irregular  acoust ic-elastic  interface  represented  by  the  surface.  S(xi). 
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Figure  3:  The  cylindrical  glass  block  model  utilized  in  this  study.  The  upper  circular  surface 
of  this  glass  block  was  etched  to  give  a  randomly  irregular  geometry  and  then  this  block 
was  submerged  in  water  to  create  an  irregular  acoustic-elastic  interface.  Measurements 
show  this  block  to  have  a  P  wave  and  S  wave  velocity  of  5600  m/s  and  3200  m/s, 
respectively.  The  density  of  the  block  is  approximately  2.65. 
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GENERATING  THE 
INTERFACE 


photoresist  layer 


Figure  4:  The  general  approach  used  to  etch  a  specific  surface  geometry  given  a  smooth 
glass  surface.  The  valley  shown  above  is  etched  in  discrete  levels  using  high  velocity 
sand  particles  and  a  photoresist  layer  which  shields  the  glass  covered  by  photoresist  from 
chipping.  Therefore  the  the  valley  is  slowly  deepened  and  broadened  at  each  level  until 
the  desired  valley  is  achieved. 
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Figure  5:  The  six  templates  used  to  generate  the  random  surface  used  in  these  experiments. 
Each  circle  has  a  19.5  cm  diameter  to  match  the  glass  surface.  Each  template  corresponds 
to  one  standard  deviation  of  depth  and  each  template  was  exposed  to  high  velocity 
particles  for  the  same  amount  of  time.  The  numbering  shows  the  order  in  which  the 
templates  were  applied. 
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(a) 


Source  Radiation  Pattern 


gure  6:  Source  radiation  pattern  for  the  two  flat-bottomed  panametrics  transducers  used 
in  this  study.  Shown  are  the  radiation  pattern  (a)  for  a  12.7  mm  (1/2  in)  source  operating 
at  1.5  MHz  (A  =  1.0  mm)  in  water  and  the  radiaion  pattern  (b)  for  a  25.4  mm  (1  in) 
source  operating  at  0.5  MHz  (A  —  3.0  mm)  in  water. 
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TQP  VIEW 


igure  7:  The  geometry  used  to  experimentally  measure  the  scattering  properties  of  a  given 
random  surface.  The  source  is  held  stationary  at  one  incident  angle  while  the  detector 
is  stepped  in  a  semi-circle  about  the  random  surface.  This  then  gives  one  realization  of 
that  surface. 


Shift  block  horizontally  by  1.25  cm. 


plane  of  receiver  1 80° 

array  (stationary) 


Shift  block  vertically  by  1.25  cm. 


'  inure  S;  I  he  rotation  scheme  used  to  generate  many  realizations  of  the  interface.  The 
sample  is  shifted  left  and  right  by  1.25  cm  steps  and  the  block  is  rotated  by  staggered 
(i()°  steps  to  give  (1  independent  realizations  of  the  surface  where  the  axes  of  rotation  for 
the  realization  sit  along  concentric  circles  on  the  surface.  This  gives  a  total  of  72  surface 
realizations.  An  additional  72  realizations  can  also  be  sampled  by  shifting  the  block  up 
and  down  by  1.25  cm  <tep>. 
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Sure  9:  Histogram  plotting  surface  height  data.  The  target  surface  height  distribution 
(a)  is  Gaussian  with  a  standard  deviation  1  mm.  The  surface  height  distribution  (b) 
based  on  profilometer  measurements  (squares)  of  the  glass  surface  is  show  i  along  with  a 
best  fitting  Gaussian  curve  (solid  line)  which  has  a  standard  deviation  of  0.6  mm.  This 
histogram  was  plotted  using  10000  surface  profilometer  measurements. 
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igure  10:  The  interface  autocorrelation  function.  The  target  autocorrelation  function  (a) 
is  a  Gaussian  function  with  a  correlation  length  e-1  of  1.0  mm.  The  actual  autocor¬ 
relation  function  of  the  glass  block  (b)  as  calculated  from  profilometer  measurements 
has  a  correlation  length  of  1.4  mm.  The  surface  profilometer  measurements  (solid  line) 
are  compared  with  Gaussian  (crosses)  and  Exponential  (circles)  functions  with  the  same 
correlation  lengths  of  1.1  mm. 
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lire  1 1:  The  surface  height  distribution  based  on  profilometer  measurements,  (a)  gives  a 
grayscale  plot  of  the  surface,  (b)  gives  a  profile  across  the  surface  as  marked  in  (a),  and 
( c )  shows  a  numerically  generated  surface  with  the  statistics  given  in  Figures  9  and  10. 
Both  (laussian  (solid  linei  and  exponential  (dashed  line)  correlation  functions  are  shown. 
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Figure  12:  Experimental  seismogram  recorded  o  the  glass  model,  (a)  shows  the  seismic 
data  recorded  as  a  function  of  angle  over  the  .  _-gular  glass  surface  with  A  =  .71a  given 
an  acoustic  beam  incident  at  20°.  The  arrow  shows  the  retroreflective  direction.  The 
source  wavelet  (b)  and  the  specular  reflection  (c)  recorded  over  a  plane  interface  are  also 
shown.  In  the  plane  layer  case  the  only  observable  energy  is  in  the  specular  direction. 
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Figure  13:  The  amplitude  recorded  experimentally  for  one  continuous  wave  realization  of 
the  fluid-glass  interface  given  an  acoustic  beam  incident  at  0°,  20°,  30°,  60°,  respectively. 
The  incident  wavelength  corresponds  to  A  =  ,71a  and  0,  is  the  angle  of  forward  (positive) 
and  back  f negative)  scattering'. 
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Figure  14:  The  2-D  mean  reflection  coefficient  (a)  calculated  numerically  over  the  Gaussian 
(solid  line)  and  exponential  (dashed  line)  surfaces  given  a  normally  incident  source  beam 
with  A  =  .Tin.  The  3-D  reflection  coefficient  (b)  recorded  over  the  experimental  interface 
(circles)  and  the  standard  deviation  of  the  finite  average  (plus)  are  also  shown.  The 
experimental  results  correspond  to  30  surface  realizations. 
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Figure  15:  Similar  to  Figure  14  except  the  incident  angle  is  now  20°  and  results  correspond 
to  20  realizations. 
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Figure  16:  Similar  to  Figure  14  except  the  incident  angle  is  now  30°  and  results  correspond 
to  10  realizations. 
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Figure  17:  Similar  to  Figure  1-1  except  the  incident  angle  is  now  60°  and  results  correspond 


to  10  realizations. 
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Figure  18:  Similar  to  Finure  12.  except  A  -  2.14«  and  results  correspond  to  10  realizations 


Figure  1!):  The  amplitude  recorded  experimentally  for  one  continuous  wave  realization  of 
the  fluid-glass  interface  given  an  acoustic  beam  incident  at  0°,  20°,  30°,  60°.  respectively. 
The  incident  wavelength  corresponds  to  A  =  2.14a  and  0S  is  again  the  angle  of  forward 
(positive)  and  back  (negative)  scattering. 


Figure  20:  The  2-D  reflection  coefficient  (a)  calculated  numerically  over  the  Gaussian  (solid 
line)  and  exponential  (dashed  line)  surfaces  given  a  normally  incident  source  beam  with 
A  =  2.1  lu.  0a  is  the  angle  of  forward  (positive)  and  back  (negative)  scattering.  The 
•5-D  reflection  coefficient  (b)  recorded  over  the  experimental  interface  (circles)  and  the 
standard  deviation  of  the  finite  average  (plus)  are  also  shown.  The  experimental  results 
correspond  to  .'50  surface  realizations. 
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Figure  21:  Similar  to  Figure  20  except  the  incident  angle  is  now  20°  and  results  correspond 
to  20  real  at  ions. 
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Figure  22:  Similar  to  Figure  20  except  the  incident  angle  is  now  30°  and  results  correspond 
to  10  realizations. 
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Figure  23:  Similar  to  Figure  20  except  the  incident  angle  is  now  00°  and  results  correspond 
to  10  realizations. 


93 


Figure  24:  The  data  recorded  about  the  retroreflection  direction  over  the  glass  interface 
using  a  piezo-film  receiver.  Shown  are  (a)  the  incident  source  amplitude  recorded  with 
the  piezo-film,  (d)  the  average  reflection  coefficient  for  65  surface  realizations  where  (b) 
the  raw  data  recorded  has  been  corrected  using  (c)  the  transmission  coefficient  of  the 
piezo-film  receiver.  The  final  reflection  coefficient  can  be  compared  to  (e)  the  numerically 
derived  reflection  coefficient  for  a  2-D  interface. 
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Figure  B.2:  The  piezo-film  receiver  with  4  layers  of  piezo-film  stacked  together  with  thin 
epoxy  layers,  (a)  The  piezo-film  sandwiched  between  two  thinner  layers  of  conducting 
glue  which  carry  current  from  the  film,  (b)  The  plane-view  geometry  of  the  receiver 
which  has  a  circular  geometry  with  a  1.25  cm  diameter,  (c)  The  parallel  connection  of 
these  piezo-elements  which  results  in  a  four  time  increase  in  current  output. 
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Figure  B.3:  Geometry  of  the  piezo-film  sheet.  The  piezo-film  element  is  inserted  into  of  a 
larger  sheet  of  stacked  plastic  layers  to  reduce  diffractions  from  the  edge  of  the  receiver 
element.  The  sheet  of  plastic  layers  w as  created  by  stacking  110 thick  layers  until  the 
impedance  contrast  matched  that  of  the  piezo-film. 
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Figure  B.4:  The  sensitivity  pattern  of  the  four  layer  piezo-film  receiver  at  500  kHz. 
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